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Preface 


5 December  1974 


V' 


This  report  is  an  analysis  of  several  signal  generation 
algorithms  which  were  proposed  and  in  some  cases  used  in  the 
Acoustic  Signal  Generation  System  (ASGS),  a multi-minicomputer 
simulation  facility  developed  at  NAVSURFWPNCEN.  In  addition 
to  deriving  the  spectrum  of  the  output  for  each  of  the  algorithms 
and  comparing  their  performance,  practical  aspects  of  implementation 
for  each  are  discussed.  The  report  will  be  of  interest  to  those 
concerned  with  the  design  of  digital  signal  generation  systems. 

The  research  reported  herein  was  partially  funded  under 
the  Digital  Acoustic  Simulator  System,  Task  Number  A370-370K/W2144-170. 

In  its  original  form,  this  paper  was  submitted  to  the  Graduate 
Faculty  of  the  University  of  Maryland  in  partial  fulfillment  of 
the  requirements  for  the  degree  of  Doctor  of  Philosophy  in 
Electrical  Engineering. 
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CHAPTER  I 
INTRODUCTION 

The  synthesis  of  specific  analog  signals  has  been 
accomplished  by  analog  systems  for  some  time.  Difficulty 
in  control  over  such  systems,  the  inherent  lack  of 
uniformity  and  reproducibility  in  their  operation,  and 
other  factors  such  as  reliability  problems  have,  however, 
limited  the  degree  of  sophistication  which  can  be 
achieved  in  analog  signal  synthesis.  In  addition,  there 
.is  a lack  of  multiplexing  capability,  so  that  many  parallel 
replicas  of  a circuit  are  generally  required.  Hybrid 
systems  employing  parametrically  programmable  analog 
devices  controlled  by  a computer  are  able  to  surmount 
some  of  the  difficulties  of  purely  analog  systems,  but 
still  have  problems  in  reliability  and  reproducibility 
of  results. 

Although  signal  synthesis  by  an  all  digital  system 
has  been  possible  for  a decade  or  more,  the  raw  computing 
power  necessary  to  perform  a reasonably  sophisticated 
Job  in  real  time  has  only  recently  become  economically 
practical.  Non  real  time  techniques  are  more  often  than 
not  impractical  from  a user’s  time  standpoint  (only 
in  rare  cases  is  one  willing  to  spend  days  generating 
minutes  of  signal  or  less),  or  completely  inapplicable 
for  on  line  requirements.  The  recent  availability  of 
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Inexpensive,  very  fast  and  very  reliable  digital 
hardware,  together  with  recent  advances  In  digital  signal 
processing  techniques  such  as  the  FFT  algorithm,  have 
resulted  in  the  feasibility  of  very  sophisticated  signal 
synthesis  by  all  digital  systems.  By  utilizing  pseudo- 
random sequence  generators  to  create  the  "random" 
components  of  signals,  completely  controllable  and 
reproducible  results  may  be  had. 

As  might  be  expected,  the  increased  versatility  and 
reliability  that  is  achievable  with  an  all  digital 
system  is  not  completely  free  from  encumbrances.  Since 
'all  representations  of  time  functions  and/or  of 
spectral  functions  must  be  discrete  and  of  finite 
duration,  aliasing  problems  are  inherent  in  every 
operation.  In  addition,  various  "clever"  implementation 
algorithms  that  appear  on  the  surface  to  be  economically 
desirable  usually  compound  the  sampled  data  system 
problems,  and  leave  considerable  doubt  as  to  whether 
the  resultant  synthesized  signal  has  the  desired 
characteristics.  It  is  these  latter  considerations  to 
which  this  paper  is  addressed. 

In  general,  signals  may  have  both  transient  components 
(finite  energy)  and  steady  state  components  (finite  power). 
This  paper  is  concerned  with  the  synthesis  of  steady 
state  components,  and  in  particular  those  that  are 
either  composed  of  one  or  more  discrete  frequency  elements. 
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or  are  random  processes  characterized  by  a power  density 
spectrum  which  varies  relatively  slowly  in  frequency. 
These  two  categories  will  be  referred  to  as  narrowband 
and  broadband  components,  respectively.  It  will  be 
assumed  that  during  synthesis  of  narrowband  components, 
the  specified  frequency  may  be  randomly  perturbed  by 
a small  percentage  to  generate  components  with  narrow 
but  non-zero  bandwidths. 

The  signals  to  be  synthesized  are  divided  into  three 
groups.  Chapter  2 discusses  the  generation  of  the 
broadband  signals  characterized  by  a finite,  reasonably 
well  behaved  power  spectral  density.  Chapter  3 considers 
the  generation  of  signals  having  discrete  or  quasi- 
discrete spectra,  and  in  particular  the  generation  of 
sums  of  such  signals  by  inverse  PFT  techniques.  As  an 
alternative.  Chapter  presents  a time  domain  technique 
for  sums  of  single  frequency  signals  whose  frequencies 
are  harmonically  related;  that  is,  all  signal  frequencies 
in  a set  are  integer  multiples  of  some  fundamental 
frequency  fQ.  In  all  cases  the  objective  is  to  develop 
and  evaluate  the  performance  of  practical,  economically 
implementable  algorithms  rather  than  to  concentrate  on 
determining  algorithms  that  are  optimum  from  a purely 
mathematical  viewpoint . 

The  glossary  on  page  vi  summarizes  the  symbols 
used  in  the  thesis <,  In  addition,  use  is  made  of  the 
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"sine"  function,  defined  by 


sine  u 


sin  nu 

tu 


since  the  form  (siniru)/iru  is  encountered  very  frequently 
in  "blocked”  data  (i.e.,  FFT)  processing.  A similar 
form  also  encountered  quite  often  is  (siniru)/(sin^u) , 
which  will  be  designated  ir,  this  paper  as  the  "sind"  function; 


sind^Cu)  * 


1 

' if 


slmru 

singu 


1.2 


which  is  3hown  in  Appendix  A to  result  from  a finite 
sum  of  exponentials 


n*0 


N*e 


4 /N-l\„ 


sindjjCu) 


1.3 


Appendix  B deriv'es  a useful  pair  of  relations  between  the 
sind  function  and  Infinite  sums  of  regularly  spaced  sine 
functions. 

It  will  be  assumed  at  the  outset  that  the  functions 
dealt  with  in  the  text  are  sufficiently  well-behaved 
to  permit  the  interchange  of  summation  and  integration 
operations  as  required,  and  other  minor  mathematical 
liberties  which  are  taken  on  occasion. 
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CHAPTER  II 

BROADBAND  SIGNAL  SYNTHESIS 
A.  General  Considerations 

The  problem  considered  in  this  chapter  is  the 
generation  or  synthesis  by  a digital  system  of  random 
signals  having  a power  spectral  density  PSD/  that  is 
controllable  over  a band  of  interest  W from  -B  to 
B Hertz.-  Since  a completely  arbitrary  control 
•capability  would  imply  an  infinite  number  of  parameters 
to  specify,  it  will  be  assumed  that  the  desired  degree 
of  detail  will  be  commensurate  with  a set  of  N or  less 
control  parameters  Yk,  k « 0,  N-l.  The  objective 

•is  the  determination  of  a signal  generation  algorithm 
defined  on  the  Yk  such  that  the  PSD  of  the  generated 
signal  y(t),  over  the  required  band  is  a reasonably 
well  behaved  function  having  the  desired  overall  form. 
As  an  example,  the  signal  generation  algorithm  could 
be  passing  white  Gaussian  noise  through  a general 
digital  filter  of  order  N,  with  the  Y^  determining 
the  pole/zoro  locations  of  the  filter. 

A further  assumption  is  made  that  the  desired  first 
order  statistics  of  the  output  signal  be  Gaussian.  This 
is  a normal  requirement  for  most  applications,  since 
natural  processes  tend  to  exhibit  Gaussian  statistics. 
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For  this  reason  It  will  also  be  assumed  that  the  digital 
synthesizer  system  hardware  Is  capable  of  generating 
samples  from  a Gaussian  process  efficiently.  Actually, 
for  most  practical  applications,  an  approximation  to 
the  Gaussian  is  adequate,  such  as  the  sum  of  4 or  more 
independent  samples  from  a uniform  process. 

If  white  Gaussian  noise  (WGN)  is  all  that  is  required 
for  the  output  signal  (power  density  spectrum  equals 
a constant  across  the  band  of  interest),  then  the  most 
straightforward  and  efficient  method  of  generating 
the  desired  signal  is  to  output  samples  directly  from 
'the  Gaussian  generator,  with  at  most  a multiplication 
by  a scale  factor  to  adjust  the  output  power  level. 

Since  the  Fourier  transform  of  a white  Gaussian  process 
is  also  a white  Gaussian  process  (exhibiting  the  proper 
symmetries  if  the  input  is  real),  specification  of  the 
signal  in  the  frequency  domain  (to  be  followed  by  an 
inverse  Fourier  transform)  requires  the  same  number  of 
independent  Gaussian  samples  as  in  the  time  domain.  Thus 
the  computational  effort  required  to  perform  the  inverse 
Fourier  transform  is  superfluous. 

If,  on  the  other  hand,  a Gaussian  signal  with  a 
colored  spectrum  is  desired,  both  time  domain  techniques 
utilizing  recursive  or  convolutional  filters  and  frequency 
domain  techniques  utilizing  the  discrete  Fourier  transform 
have  their  application.  If  a relatively  simple  spectral 
shaping  is  required,  passing  the  white  Gaussian  (time) 
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samples  obtained  directly  from  the  generator  through  a 
time  domain  filter  may  be  adequate.  This  approach,  however, 
has  two  rather  severe  limitations.  The  computational  load 
exceeds  the  equivalent  FFT  implementation  load  if  the  filter 
function  required  has  more  than  about  5 to  10  poles  or 
zeroes.  Worse  yet,  design  of  a time  domain  filter  to 
implement  a given  frequency  domain  specification  is  not 
straightforward,  and  may  become  impractical  if  the 
spectral  shaping  must  be  varied  in  real  time.  The  time- 
domain  filtering  approach  is  practical  then  only  if  very 
simple  shaping  is  desired,  or  most  profitably,  if  the 
desired  shaping  is  naturally  specified  in  terms  of  a few 
recursive  filters  (high,  low  or  bandpass  functions  with 
rolloff  in  multiples  of  6 db/oct.)  Computational  load 
for  convolutional  filters  becomes  excessive  beyond  a few 
zeroes,  so  that  the  required  convolution  is  often  per- 
formed by  transforming  to  the  frequency  domain, 
multiplying,  and, transforming  back.  This  of  course 
becomes  the  same  thing  as  frequency  domain  filtering. 

If  considerable  detail  in  the  spectral  power  density  is 
required,  frequency  domain  filtering  techniques  utilizing 
an  FFT  processor  become  quite  attractive.  The  desired  power 
spectral  density  can  be  specified  at  N equally  spaced  points 
from  -B  to  B in  the  frequency  domain,  or  if  an  impulse 
response  or  autocorrelation  function  is  given,  the  equivalent 
spectral  weighting  function  can  be  easily  obtained  by 
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a forward  transform^  Once  some  threshold  value  of  N 
has  been  passed  and  the  decision  to  use  FFT*  techniques 
has  been  made,  the  computational  load  is  only  weakly 
dependent  on  degree  of  detail  implemented.  The  major 
disadvantage  of  the  FFT  methods,  though,  is  the  fact 
that  since  processing  must  be  done  in  blocks  rather 
large  amounts  of  memory  are  required.  The  memory 
requirements  are,  in  fact,  proportional  to  the  degree 
of  detail  required,  i.e.,  the  parameter  ii, 

B.  Frequency  Domain  Filtering  Algorithm 

The  most  straightforward  approach  to  synthesis 
of  broadband  signals  by  FFT  is  to  implement  directly 
the  frequency  filtering  diagram  shown  in  Figure  2.1. 

First  white  Gaussian  noise  is  generated  in  the  time 
domain  as  before.  The  forward  transform  of  a block 
of  N samples  is  computed,  yielding  N/2  single-sided 
complex  frequency  coefficients.  These  are  then 
multiplied  by  the  corresponding  coefficients  of  the 
filter  weighting  function.  The  inverse  FFT  of  the 
result  yields  the  desired  filtered  time  samples. 

However,  to  avoid  erroneous  results  due  to  circular 
convolution,  the  actual  filtering  implementation  must 
incorporate  either  of  the  ,:overlap-save"  or  the  "overlap- 
add”  computational  techniques. ^ In  either  of  these 
schemes,  the  desired  filter  Impulse  response  of, 
for  example,  N non-zero  samples  must  be  augmented  by 
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a block  of  K zeroes  to  produce  a total  block  size  of 
M * N + K time  samples.  The  forward  transform  of  this 
function  produces  the  spectral  filter  function  desired 
having  M/2  single-sided  complex  coefficients.  Input 
data  in  blocks  of  M samples  (or  N samples  plus  K zeroes 
for  the  overlap-add  version)  are  then,  forward  transformed 
to  produce  M/2  spectral  values.  After  multiplication 
of  the  input  spectra  block  by  the  augmented  filter 
weighting  function  and  inverse  transforming,  K samples 
of  the  output  are  free  of  circular  convolution  error. 

The  price  to  be  paid,  however,  is  increased  memory 
size  by  K words,  and  a required  overlapping  of  the  input 
blocks  by  N samples. 

The  complete  process  is  illustrated  in  Figure  2.2 
for  the  overlap-save  method.  Note  in  2.2b  that  to 
specify  N parameters  for  the  spectral  detail  (N/2 
complex  coefficients)  a block  size  M » N + K must  be 
used  in  performing  the  transforms.  For  computational 
efficiency,  it  would  intuitively  seem  from  2.2a  that 
K should  be  considerably  larger  than  N.  However, 
memory  requirements  are  directly  increased  by  increasing 
K,  and  from  the  memory  standpoint  K should  be  made  as 
small  as  possible.  In  addition,  the  computational 
gains  are  almost  non-existent  for  K much  greater 

u 

than  N,  (due  to  the  log  factor  in  the  (j-log  M)  FFT 
load  formula)  and  actually  become  losses  for  M more 
than  about  8N.  Figure  2.3  illustrates  the  computational 


10 


NOLTR  74-215 


M = N+K  SAMPLES  FOR  FIRST  TRANSFORM 


M = N+K  FOR  SECOND  TRANSFORM, ETC 


2.2a  INPUT  WAVEFORM  BLOCKING,  N SAMPLE  OVERLAP 


H(K) 

N/2 

COMPLEX 


uJ 

\ 

IDFT 

Drr 

^ H'(K) 

N K ZEROES  ADDED  ^2  =(N  + K)  12 

REAL  COMPLEX 


2.2b  DERIVATION  OF  REDUNDANT  FILTER  COEFFICIENTS 


FILTER  COEFFICIENTS 
M/2  COMPLEX 


INPUT 
BLOCKS 
M REAL 


INPUT 
SPECTRA 
M/2  COMPLEX 


FILTERED  OUTPUT 
M REAL 

K VALID  SAMPLES 
N SAMPLES  DISCARDED 


2.2c  OVERLAP-SAVE  BLOCK  DIAGRAM 


FIG.  2.2 


OVERLAP-SAVE  FILTER  ALGORITHM 
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MEMORY, 

LOAD 


N = 32 


N = 512 


£ RELATIVE  MEMORY  SIZE 
X RELATIVE  COMPUTATIONAL  LOAD 


FIG.  2.3  COMPARISON  OF  MEMORY  REQUIRED  AND  COMPUTATIONAL  LOAD 
VS  M FOR  OVERLAP-SAVE  ALGORITHM 
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load  and  memory  required  as  a function  of  M for  N * 32 
and  N * 512.  For  each  plot,  "i"  is  the  requirement 
for  the  base  transform  (32,  ,512).  In  both  cases  M * 2N 
appears  to  be  a good  compromise  between  FFT  load  and 
memory  required. 

Assume  that  the  power  spectral  density  (PSD)  $(f)  has 
been  specified  at  N-l  equally  spaced  points  covering 
the  band  of  interest  W,  i.e., 

#(kfc)  - Yk  2;  k » - | + 1,  ...»  0,  ...,  N/2-1;  2.1 

where  N is  assumed  even,  Yk*Y_k  are  real  and  positive,  and 

where  f , the  resolution,  is  given  by 
c 

f„  « W/N  2.2 

v 

The  Yk  are  then  the  magnitudes  of  the  amplitude  response 
at  the  frequencies 

fk  “ 2.3 

The  spectral  response  in  the  band  of  interest  W 
at  frequencies  lying  between  the  fk  is  determined  by 
the  particular  synthesis  algorithm  employed  to  generate 
the  output  time  sequence  yn.  For  the  overlap-save  filter 
algorithm  described  above,  assume  that  the  Yk  are  used 
directly  as  the  N-l  non-zero  spectral  coefficients  of 
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the  filter  function,  i.e.. 


Hk  " Yk»  k “ "N/2  + 1»  N/2  - 1 2.4: 


where  H_k  « H#^  * Y^  since  the  Yk  are  real.  The  correspond 
ing  time  sequence  is  given  by 


n 


N/2-1 

I , Hk  e 
k«-N/2  K 


j2imk 


N/2-1 

l Y.e  w ; n=0, . . . ,N-1 
k*-N/2  K 


2.5 


where  Y is  assumed  to  be  zero. 

The  filtering  operation  to  be  performed  is  then 
the  convolution  of  a function  h(t)  with  the  input  function 
x(t),  with 


m ** 

x(t)  *=  l xm  «(t  - mTg)  2.6 

- m=»-«  ' 

where  the  xm  are* the  white  Gaussian  noise  samples,  and 
where 

N/2-1 

h(t)  = l h 6 (t  - nT  ) 2.7 

n=-N/2  n s 

1 

and  Ts  = jjjr-  is  the  basic  sample  period  of  the  system. 

To  determine  the  spectrum  of  h(t),  consider  a function 
h (t)  which  is  periodic  with  period  NT  - 1/f  and  is 

r o v 

identically  equal  to  h(t)  over  the  interval  [-jTs,  jyTg], 

A function  which  consists  of  a repeating  or  periodic 
sequence  of  impulses  has  a spectrum  that  is  also  a 
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periodic  sequence  of  impulses,  where  the  Impulse  spacing 
in  each  domain  is  just  the  reciprocal  of  the  repeat 
interval  in  the  other  domain.  The  Impulse  values  over 
one  repeat  interval  in  each  of  the  two  domains  are 
related  by  the  Discrete  Fourier  Transform  (DFT) 
modified  by  a scaling  factor  of  T.  Since  by  definition 
the  hn  are  the  inverse  DFT  of  the  H^,  the  spectrum  of 
hp(t)  is  given  by 


m 

i 

m*-» 


N/2-1 

k«-N/2 


fljr  Hk  «(f  - (k+mN)fc) 
8 


Thus  the  spectrum  of  hp(t)  is  sampled  and  periodic  with 
the  being  proportional  by  1/NT8  to  the  sample 

values  over  one  period  Nf_.  But  h(t)  is  just  h (t) 

c p 

multiplied  by  a unit  pulse  of  width  NT.  • r , 

c 

Therefore  the  spectrum  of  h(t)  is  Hp(f)  convolved  with 
jr-  sine  f/fc,  or 

— Hk  sincC  (f  - (k+mN)f  )/f  ]' 
NTS  K c c 

2.< 


\ H(f ) 


1 

77 


c m*-» 


N/2-1 

k«-N/2 


Interchanging  sums  and  rewriting  slightly, 

NT  N/2-1  - f 

H(f)  * twr”  I I sincC*—  - k - mN]  2.] 

WJ,s  k«-N/2  K m»-»  *c 


or,  using  B19, 


f- ' 

-1 
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N/2-1  r 

H(f)  - l Hk  sindN(J-  - k)  2.11 

k«-N/2  k N rc 

j 

i Since  the  input  x(t)  is  assumed  to  be  white, 

«*> 

| its  PSD  *x(?)  is  just  some  constant  which  can  set  to 

i *,  . 

| unity  for  convenience.  The  PSD  of  the  signal  y (;t) 

i 

j at  the  output  of  the  filter  is  then 

♦ M)  * ir(f)  • |H(f ) |2  2.12 

; yA 

or 

N/2-1  f 2 

♦ ,(f  ) * l Hk  sindM(~-  - k)  2.13 

yx  k*-N/2  K N rc 

Over  the  positive  spectrum  from  0 to  B Herts,  is 
Just  equal  to  Y^,  and  the  spectrum  of  the  output  over 

f 

the  band  of  Interest  is  given  by  2.13.  Thus  the  expected 

! value  of  the  spectral  magnitude  of  the  signal  generated 

i * 

1 by  the  above  convolution  filtering  algorithm  passes 

through  the  values  Y^  at  the  frequencies  kfc,  k « -N/2, 

...,  N/2  - 1 since  sindN(u)  * 0 for  non-zero  integral 
values  of  u < N and  exhibits  sampling  function 
interpolation  at  points  between  the  kfc.  Note,  however, 
that  the  interpolation  includes  contributions  from  the 
replicas  of  the  outside  the  band  of  interest,  and  thus 
results  in  the  interpolation  form  sind^Cu). 
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C.  Concatenated'  Segment  Algorithm 

It  would  seem  that,  since  the  transform  of  a V6N 
time  function  is  a WGN  spectrum  (reals  and  lmaginarys 
in  the  spectrum  are  all  independent  Gaussian  random 
variables),  it  would  be  possible  to  eliminate  the  input 
forward  transform  and  to  generate  the  WGN  spectrum 
directly.  The  problem  is  that,  although  within  each 
block  the  N random  variables  are  Independent  and 
Gaussian  in  both  domains,  there  is  a required  correlation 
between  blocks  due  to  the  N overlapped  samples. 

It  would  be  necessary  to  generate  the  spectral  random 
variables  such  that  the  last  N samples  of  the  transform 
f of  one  are  identical  to  the  first  N samples  of  the  transform 

^ of  the  next.  It  Is  not  immediately  obvious  how  to 

generate  spectra  with  the  required  constraint,  since 
the  time  domain  is  the  "natural”  domain  for  specifying 
the  phenomena  desired. 

j If  circular  convolution  is  completely  ignored,  the 

t computational  load  can  be  reduced  from  one  forward  plus 

one  inverse  transform  per  output  block  multiplied  by 
some  factor  for  overlap,  to  Just  an  inverse  transform 
; with  no  overlap.  Assuming  an  overlap  of  two  for  the 

fast  convolution  algorithm,  this  corresponds  to  a 
potential  reduction  factor  of  four  in  computational 
loading.  It  is  thus  tempting  to  try  to  obtain  satis- 
I factory  results  by  implementing  the  following  signal 

: generation  algorithm. 
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Let  Hk,  k « 0,  1,  . ..,  N/2  - 1 be  the  single-sided 
spectral  coefficients  of  a discrete  filter  weighting 
function.  Let  a WON  spectrum  be  generated  with  N/2 
complex,  single-sided  frequency  coefficients,  X^,  where 
the  probability  density  of  both  the  real  part  x£  and 
the  imaginary  part  x£  of  X^  is  given  by 


2 

p r(X)  * p ±(X)  - -L-  e“X 

Xk  ?k 


2.14 


2 1 

i.e.,  a zero-mean  Gaussian  process  of  variance  a B p. 

r*  4 

The  N random  variables  x£,  Xj^j  k = 0,  1 ...  N/2,  are 
all  independent.  The  magnitudes  of  the  X^  are  then 
Rayleigh  distributed 


PJX  | ^ = 2re"r  > r L 0 2.15 

. k 


- ~~2 
with  mean  r * /*  /2  and  mean  square  r =1, 

The  phase  angles  of  X^  are  uniformly  distributed 

over  2n  radians: 


p x (® ) * 0^9  < 2n  2.16 


The  signal  generation  algorithm  consists  of  first 
forming  the  single-sided  spectrum 

Yk  s VXk;  k B °*  *.»  N/2  - 1.  2.17 
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Then  an  inverse  Discrete  Fourier  Transform  (XDFT)  is 

performed  on- the  Yv  to  produce  the  real  time  series  ,.y  : 

n * 0,  ...  N - 1.  The  yn  become  the  sample  values  of  -the 

output  broadband  process  y (t)  for  an  interval  NT  . where 

T_  is  the  basic  system  sample  period  as  before.  The 
s 

algorithm  is  repeated  indefinitely  with  a new.  Independent 

f*  Vi 

set  of  X^’s  on  each  iteration.  The  l iteration  generates 
the  set  of  X’s,  X^,  from  which  the  samples  of  yx(t) 

for  the  interval  (1-1/2 )NT_  < t < (£+l/2)NT_  are  obtained. 

S **  5 

Thus: 

v « N/2-1 

yX(t)  * l l y > 6(t  [£N+n]T  ) 2.18 

n=-N/2  s 

where  yn^  is  the  value  of  yn  generated  on  the  l n iteration. 

Each  segment  of  the  output  function  generated  in  a 

single  iteration  can  be  thought  of  as  exactly  one  period 

of  a periodic,  sampled  function  yp£(t),  with  sample  period 

T_  and  fundamental  period  NT_.  The  spectrum  of  such  a 

function  is  also  a sampled,  periodic  function  h i(t)- 

with  period  f = <1/T_  and  sample  period  f = 1/NT  . 

Since  y p^(t)  and  Yp^(f)  are  related  over  one  period  of 

each  by  the  Discrete  Fourier  Transform,  Yp£(mfc)  is  just 

equal  to  Y^  • .for  m = k + rNj  m and  r integers. 

s 

Thus 


Nr 


i 


-P4  k--N/2  Nis  c 


th 


2.19 


with  being  the  Y.  generated  on  the  l iteration. 

i Kv  *.  it 
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The  power  density  spectrum  * (f ) of  the  function 

yp<e(t)  is  Just* 


*v  (f)  = '—±ry  l 1 l |Yk.|2  «(f-[k+rN]f  ) 2.20 

-yp l (NT  y k=N/2  r=-~  ~K4  c 

5 


The  autocorrelation  of  the  final  output  y (t)  can  be  found 

by  averaging  over  time  the  statistical  autocorrelation 

function  <{>  (t,x).  The  latter  is  just  the  autocorrelation 

function  of  the  member  functions  yp^  if  t + t is  still 

in  the  same  segment  as  t and  zero  otherwise,  A (t) 

yp  l 

is  obtained  by  first  taking  the  ensemble  average  of 
, and  then  performing  the  Fourier  transform. 

~yp  Z 

Applying  the  expectation  operator  to  2.20,  only  the 
Yk£  are  affected  and 


N/2-1  *> 

* (f)  = — 1 l l y 2 6 (f-[k+rN]f  ) 2.21 

yp£  (NT  )2  k«-N/2  r=-»  k c 

5 


since 


■ V 


2.22 


*The  interested  reader  can  readily  convince  himself 
of  the  validity  of  2.20  by  considering  the  Power 
Spectrum  represented  by  2.19  and  noting  the  independence 
of  the  Yk£  for  all  k and  l . 


20 


NOLTR  74-215 


by  definition.  The  time  average  of  4 (t,x)  is  then  Just 

<y 

$ (t)  multiplied  by  the  percentage  of  all  possible 

yp  l 

times  t that  a delay  of  t still  remains  in  the  same 
segment.  The  multiplying  function  is  thus  triangular 
with  a value  of  unity  at  the  origin  (average  power 
in  yx(t)  is  same  as  average  power  in  member  functions) 
and  decreases  linearly  to  zero  at  + the  segment  length 
NTs . Thus 

0 < x < NT 
— s 

M > NTg  2.23 

-NTg  < t < 0 


$ (OCl-jw-]; 
yp£  Nis 


♦y(l)  “ 


*V(,Ul+*7]; 


The  output  power- density  spectrum  is  then  the 
convolution  of  * (f)  with  the  transform  of  the 

yP-e 

triangular  pulse.  The  latter  is 


PCD  = 


-X*  sine2  £- 

V 


2.24 


Convolving  with  2.21,  the  PSD  of  the  output  is 


(NT  V 

*v(f)=  S 


(NTs)' 


N/2-1  5 » 

l V*  i 

k=-N/2  K r=-« 


sinc2(X 

c 


k -rN) 
2.25 
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or,  using  B24 


N/2-1  2 2 f 

* s l V sindM‘ £tfr  - k)  2,26 

yx  k=-N/2  K N xc 

Comparing  2.26  with  2.13,  the  PSD  differ  in  that  the 
PSD  of  the  convolution  filter  algorithm  is  of  the  form 
"square  of  sum",  whereas  the  segment  function  algorithm 
yields  the  form  "sum  of  squares".  Since  the  sind 
functions  have  zeroes  at  all  the  cell  frequencies  except 
the  source  of  each,  both  algorithms  yield  a PSD  that 
passes  through  the  points  | | at  frequencies  kfQ. 

The  difference  exists  only  for  those  frequencies 
lying  between  the  kf  . The  contributions  to  these 

V 

frequencies  are  due  to  a continuous  impulse  response 
function  for  the  convolution  filter.  This  response 
function  i3  a sum  of  sind  functions,  but  each  is  excited 
by  the  same  random  proce.ss,  and  therefore  the  appearance 
of  the  "coherent  sum"  form.  On  the  other  hand,  however, 
each  segment  output  for  the  second  approach  is  obtained 
by  summing  a set  of  independent  narrowband  processes, 
each  with  a spectrum  of  the  form  sindj.(x)  and  as  is 
expected  these  contributions  add  incoherently  to  form 
the  final  output.  Hov/ever,  since  the  stated  requirement 
was  to  produce  broadband  noise  with  a PSD  that  is  well- 
behaved  and  continuous,  and  that  passes  through  the  points 
jHjJ  at  frequencies  kfQ , both  algorithms  appear  to 
qualify.  Thus  if  the  PSD  is  the  sole  criterion  to  be 
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satisfied,  a considerable  savings  In  computational 
load  can  be  realized  by  implementing  the  simpler 
algorithm. 

However,  all  is  not  well  with  the  "short-cut" 

approach.  Each  segment  of  the  output  function  is 

actually  one  period  of  a periodic  waveform,  and  as 

such  is  continuous  from  one  cycle  to  the  next,  or 

equivalently,  from  the  end  to  the  beginning.  The  next 

segment  generated  in  this  manner  is  completely 

independent  of  the  first,  and  the  value  of  the  function 

and  its  derivatives  at  the  end  points  .will  have  no 

relationship  to  the  coresponding  values  at  the  ends 

of  the  adjacent  segments.  Thus  the  penalty  for  ignoring 

the  circular  convolution  problem  is  a discontinuity 

in  the  function  and  its  derivatives  every  N output 

samples.  It  is  interesting  to  note  that  the  first 

order  density  function  of  the  output  yx(t) 

(evaluated  at  the  sample  points)  is  independent  of 

/ 

position,  i.e., 


PyxfXnTg)]  = p'x(y)  £ — — 
y /2?c2 


2.27 


where  the  "approximately  equal"  becomes  an  equality  if 
the  reals  and  imaginarys  are  truly  Gaussian. 
However,  although  the  power  spectrum  of  the  generated 
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signal  is  also  as  desired,  the  fact  that  the  time 
function  is  composed  of  statistically  Independent 
segments  is  apparent  in  the  second  order  density 
function.  Although  the  temporal  autocorrelation 
function  of  the  output 

R(t)  * / yx(t ) yx(t+t )dt  2.28 

(the  transform  of  the  PSD)  has  no  indication  of  the 
anomaly,  the  ensemble  correlation  function 

♦ (t,r)  « E{yx(t,)yx(t+x)>  2.29 

y 

is  not  independent  of  t and  therefore  is  not  equal  to 
R(x ) . Thus  the  output  process  is  non-ergodic.  The 
practical  implications  depend  on  the  application, 
but  in  general  the  simple  algorithm  would  not  likely 
be  satisfactory  if  the  output  is  to  be  fed  into  a 

phase  or  transient  sensitive  processing  system  (such  as 
the  human  ear). 

D.  Overlapped  Segment  Algorithm 

For  applications  in  which  the  discontinuities  in 
the  output  signal  are  unacceptable,  various  compromise 
algorithms  are  possible  which  require  overlapping 
output  segments  and  therefore  increased  computational 
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load,  but  which  at  least  result  in  continuity  of  the 
function.  Consider  again  the  overlap-add  version  of 
the  convolution  filter  algorithm  (see  reference  (1), 
p.  210)  with  an  overlap  of  two,  i.e.,  M ■ 2N.  The 
filter  impulse  response  and  N points  of  the  input 
waveform  are  both  augmented  by  K=N  zeroes.  The  spectra 
of  the  augmented  functions  have  N complex  coefficients 
over  the  same  band  that  the  unaugmented  functions 
have  N/2 j the  additional  N/2  coefficients  representing 
spectral  points  halfway  between  the  original  points. 

The  new  points  represent  interpolated  values  based 
on  the  sind^(x)  weighting  intrinsic  to  the  discrete 
Fourier  transform.  Note  that  this  implies  that  each 
interpolated  value  is  obtained  from  a weighted  average 
of  all  N/2  of  the  original  spectral  coefficients. 

It  is  obvious  that  performing  this  interpolation 
algorithm  in  the  frequency  domain  would  require 
considerably  more  computation  than  the  FFT  of  the 
original  input  time  function.  (The  filter  response 
only  needs  to  be  transformed  once  at  setup,  so  it 
is  not  a problem).  However,  it  is  possible  to  apply 
other  less  sophisticated  interpolation  algorithms, 
which  in  general  result  in  an  approximation  to  augmenting 
the  input  with  N zeroes.  Thus  some  improvement  over 
the  concatenated  segment  algorithm  might  be  had  by 
overlapping  and  combining  segments  that  are  derived  by 
a simple  interpolation  scheme. 
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The  above  considerations  suggest  an  algorithm 
based  on  an  overlap-add  scheme,  .but  where  the  segments 
overlapped  are  still  Independent.  Let  segments  of 
length  N be  generated  as  in  the  previous  algorithm, 
but  then  be  windowed  or  modulated  and  overlapped 
before  combining.  To  prevent  the  output  signal  from 
exhibiting  modulation  related  to  the  segment  period, 
the  nodulated  last  half  of  segment  l must  add  to  the 
modulated  first  half  of  segment  M 1 to  yield  first 
order  statistics  that  are  Independent  of  position. 

Since  the  random  process  for  each  segment  is  Independent 
of  all  other  segments,  the  overlapping  halves  add 
incoherently.  Thus  an  appropriate  modulating 
function  is  a sine  or  cosine  half  cycle  pulse  that 
has  zeroes  at  each  end  of  the  segment  and  is  unity 
at  the  center.  The  interpolation  algorithm  that 
produces  a segment  modulated  by  such  a sine  pulse  is 
easily  derived. 

Let  yp^(t)  be  the  sampled,  periodic  segment  function 
derivable  from  the  Yj^'s  by  the  Discrete  Fourier  Transform. 

A 

It  is  required  that  the  spectrum  of  S^Ct)  be  generated; 
where 

yp.e(t)  s * cos  » 2*30 
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and  where  T_  -is  the  basic  sample  rate  as  before. 
Multiplication  of  the  time  functions  corresponds  to 
convolution  of  the  transforms.  The  Fourier  Transform 
of  cos  jpjr-  is  Just 


Mtr)  - |c*cr  - 5^-0  + iff  + ji-)] 

s s 


2.31 


i.e..  Impulses  at  +1/2NTS.  The  cell  spacing  is  1/NTS, 
so  M(f ) consists  of  impulses  at  times  the  cell 
spacing.  Obviously  since  the  function  Y »(f)  (spectrum 

A ^ 

°f  yp^(t ))'  is  to  be  non-zero  for  only  multiples  of 
fc  * l/NT^,,  direct  application  of  2.31  is  not  useful. 
However,  by  applying  a little  sleight-of-hand,  it  is 
possible  to  obtain  the  spectrum  of  a signal  that  is  / 

functionally  equivalent  to  yp^(t)  and  that  has 
components  only  at  f ■ kfc  for  integer  k.  The  trick 
is  to  redefine  the  Hk’s  as  the  spectral  magnitudes 
at  frequencies  f = (k  + |0fc.  Convolution  of  each 
of  the  spectra  Y^(f)  generated  from  the  new  Hkfs 
with  M(f)  yields  the  desired  spectrum  where 

A 

each  component  ¥p_£(kfc)  is  the  average  of  the  two 
adjacent  components  of  Y^Cf),  Yk^  and  Y(k+1)^. 

A 

The  transform  of  Yp^(f)  then  produces  the  function 

A 

yp*(t ) , of  which  one  NTg  - second  segment  between 
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zeroes  of  the  sinusoidal  envelope  is  added  to  the  output 

yx(t)  being  generated.  The  fact  that  the  phases  of  the 
•* 

* 

components  of  Y^Cf)  are  uniformly  distributed  over  2* 

radians  permits  taking  the  NT  -second  segment  of  y »(t) 

-NT  NT  s 

from  to  -sp  with  no  loss  in  generality.  To  find  the 
PSD  of  the  output  generated  by  the  latter  method,  let 

Yp^(f)  be  written 

• N/2-1  ..  . , 

• I I nr-  *k,  «(r-(rN+k+|)r  ) 2.32 

- p 1 r— « k-N/2  NTs  M i 0 

where  the  are  generated  from  the  redefined  H^’s  as  in 

2,17.  Convolving  with  M(f) 


N/2-1 

YD/(f)‘2 UT-  l l Yk/{«Cf-(rN+k)f  + «(f-(rN+k+l)f  } 
-pt  *N1S  k=_N/2  c c 

2.33 

Isolating  a segment  of  length  NT  results  again  in 
convolving  with  ^ sinc  f/f  c » yielding ' 


. • N/2-1 

Y,(f)  st  I I Ykf{sinc(f/f-rN-k)+sinc(f/fc-rN-k-l)} 
-L  2r=-«  k*-N/2  c 

2.31* 

or 


!£(f> 


. N/2-1  - 

? E Yk-t  I {sine (f/fc-rN-k)  + sine (f/f c-rN-k-l ) } 
* ka-N/2  ' r=-“ 


2.35 
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* 

The  Energy  Spectral  Density  of  7^ CD  is  then 


#y/f)  * EC!f(f)  YICf)J 


2.36 


, M/2-1  M/2-1  _ f 

1 l l E[Y  .Y7J-  [ [sinc(iir — rN-k) 

2 k*-N/2J*-N/2  J r—  c 


+ sine  ( jr*  -rN-k-1 ) 2 * X C sine  ( jr-  -mN-J  )+sinc  ( jr-  -mN-j  -1 ) ] 


m»-» 


2.37 


Again,  is  independent  of  unless  k«J,  and  then 

E«w  ^ ■ 4 

■l  N/2rl  «■  * ft  f>  j 

K (O-i  I K ( l (sine (4-  -rN-k)+sinc (4-  -rN-k-1)})2 

yl  2 k*-N/2  K r»—  *c  xc 

2.38 

Finally,  the  PSD  of  yx(t)  is  again  just  twice  the  average 

«* 

over  t of  2.38  (factor  of  two  due  to  overlap),  or 


♦ x(f ) 

y 


2*#;  (f) 

*1 


2.39 
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The  form  of  2,38  is  the  same  as  2.10  and  2.25  except 
that  the  interpolating  function  is  the  sum  of  two 
sinc(x)  functions  displaced  by  unit  separation. 

Writing  the  sum  out, 

sine  * + sine  (x  + 1)  - SJOi  + ll 

' ^x-TT'  2-« 

Note  that  the  result  has  the  same  oscillatory  behavior 

2 

as  sine  x,  but  falls  off  as  x rather  than  x.  This  is 
especially  beneficial  when  it  is  desired  to  produce 
an  output  yx (t ) whose  spectrum  has  abrupt  changes  in 
level  as  a function  of  frequency. 

Summarizing,  the  straight forward  application  of 
standard  convolution  filter  “echniques  to  the  synthesis 
of  a random  process  with  specified  PSD  requires  on 
the  order  of  a factor  of  four  more  computational  load  than 
a minimum  implementation  ignoring  circular  convolution. 
Although  the  PSD  at  the  specified  points  and  the  first 
order  time  statistics  are  .the  same,  the  minimum  implementa- 
tion  output  is  discontinuous  and  may  be  unsatisfactory 
for  transient-sensitive  systems.  A compromise  solution 
requiring  a factor  of  two  overlap  and  therefore  a factor 
of  two  more  computational  load  is  somewhat  more  satisfactory 
in  that  the  output  is  continuous,  but  the  derivatives 
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at  the  segment  ends  are  still  discontinuous.  However, 
it  will  be  shown  in  Chapter  3 that  the  factor  of  two 
overlap  is  very  desirable  when  discrete  or  line 
components  are  to  be  added  to  the  generated  signal. 
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CHAPTER  III 

IPFT  GENERATION  OP  DISCRETE  FREQUENCIES 
A.  The  Exact  Synthesis  Algorithm 

In  an  analog  system,  the  generation  of  a sinusoidal 

waveform  is  usually  implemented  by  constructing  an 

oscillator  having  the  desired  frequency  . 

The  direct  adaptation  of  this  approach  to  discrete 

or  line  frequency  generation  by  a digital  system  requires 

a second  order  difference  equation  with  a pole 

a't  the  desired  frequency.  As  a general  digital  synthesis 

technique,  this  approach  has  problems  with  stability 

due  to  a coefficient  quantization,  as  well  as  the 

disadvantage  that  each  frequency  requires  a separate 

generator.  On  the  other  hand  samples  of  the  desired 

sinusoid  are  generated  one  at  a time,  and  thus  for  a 

small  number  of  discrete  frequencies  desired,  the 

# 

meager  storage  requirements  are  attractive  relative  to 
block  processing  techniques.  However,  for  the  case 
where  the  sum  of  many  discretes  is  desired,  significant 
computational  savings  can  be  obtained  by  utilizing 
block  processing  with  the  Fast  Fourier  Transform  (FFT). 

Implementation  of  a general  second-order  difference 
equation  requires  about  the  same  amount  of  computation 
as  an  FFT  butterfly.  Generation  of  N samples  of  a 
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single  sinusoid  by  difference  equation  requires  N 

such  computations,  whereas  generation  of  N samples 

by  IPFT  requires  ^log^  computations  for  the 

transform  plus  whatever  computation  is  required  to 

specify  the  desired  sinusoid  in  the  DPT  domain. 

Thus,  for  N on  the  order  of  1000,  generation  of  the 

sum  of  more  .'than  5 discretes  begins  to  favor  the 

FFT  approach  unless  an  unreasonable  amount  of 

computation  is  required  to  specify  the  lines  in  the 

sampled  frequency  space  input  to  the  IFFT. 

For  generation  of  discretes  with  frequencies  that 

are  exactly  equal  to  the  cell  frequencies  fc  * fg/N, 

where  f is  the  sampling  frequency,  the  specification 

in  the  sampled  frequency  domain  is  trivial.  Consider 

the  IFFT  of  N/2  single-sided  complex  frequency  coefficients 

Fk,  k=0,l, . . .N/2-1,  to  produce  N time -samples  f , 

n=0,...N-l.  A non-zero  Fk,  say  , will  result  in  a 

contribution  to  the  output  time  segment  of  exactly 

cycles  of  a sinusoid  (period  = N/^  sample  periods), 

with  a magnitude  and  a phase  relative  to  the  beginning 

of  the  segment  determined  by  the  complex  coefficient 

Ft  . Since  successive  segments  generated  with  the  same 
K1 

F.  .fill  match  perfectly  at  the  segment  boundaries, 

the  result  will  be  a perfect  sinusoid  (sampled,  of 

^1 

course)  with  frequency  -yj-f g , Thus  for  independent 
discretes  with  frequencies  in  multiples  of  fg/N,  a 
single  complex  number  added  to  the  appropriate 
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spectral  cell  before  each  inverse  transform  will 
result  in  the  desired  sinusoid  added  to  the  output 
time  waveform. 

For  frequencies  that  are  not  exact  multiples  of 
f_/N,  the  situation  is  considerably  less  straight- 

D 

forward.  Since  the  discrete  Fourier  Transform  (DFT) 
is  information  lossless  and  reversible,  the  easiest 
way  to  see  what  is  required  to  generate  arbitrary 
frequency  sinusoids  is  to  consider  the  forward 
transform  of  successive  segments  of  the  desired 
waveform.  First  consider  the  (complex  to  complex)  DFT 
of  an  arbitrary  complex  exponential  w(t)  given  by: 

j (2uf  t+a) 

w (t ) = e a 3.1 

where  f is  between  0 and  f/2,  and  o is  the  phase  at  the 

d 5 

beginning  of  the  segment  of  length  NT_  'from  which  N 

s 

samples  have  been  taken.  Let  f_  be  expressed  as 

a 

' fs 

fa  - (m  + d)-f  3.2 

fs 

where  -w-  is  the  cell  spacing,  and  m is  an  integer  and 

• fs 

d less  then  unity,  m-jj-  is  then  the  next  lowest  cell 

frequency,  and  d is  the  fraction  of  a cell  spacing 

■ fs 

that  f is  above  m-jp  . The  sampled  version  of  w(t) 
is  then 

« ■ K 

W*<t)  = l eJaeJ2',Cm+d)-r  nTs  s(‘-nTs)  3.3 

n--“ 
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and  the  samples  of  the  segment  of  interest  are 

J (njrCm+d  )n+o ) 

w(n)  ® e , n^O,  1 ...  N-l 


3.4 


where  the  fact  that  f T « 1 has  been  noted.  Thje  DPT 

s s 

is  then 


W(k) 


1 V J.  fJT",d)  -JT 

J I . e 6 e 


n«0 


3.5 


Removing  the  phase  angle  from  the  sum  and  combining 
the  exponential  arguments 

.Sc,  N-l  -J^(k-m-d) 

W(k)  ' V I e N 

n=0 


3.6 


From  Appendix  A, 


N-l  -jajnu 

l e w = Ne  sindM(u) 

n=0 


3.7 


and 


W(k)  «=  eJ  e 


ja  -jJCN-l)(k-m-d) 


sindN(k-m-d) 


3.8 


If  d=0,  only  W(m)  » e^a  is  non-zero  as  was  noted  above. 
However,  for  any  non-zero  d,  0 < d < 1,  all  W(k)  are 
non-zero.  The  magnitude  of  the  W(k)  are  governed  by 
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the  sindM(k-m-d)  form,  and  the  phases  are  governed 
by  e+^  “ * 5 Ck-m-d) ) ^ where  0 is  the  desired  phase 

angle  at  the  beginning  of  the  segment. 

For  a real  tim</ function 

w(t)  + •>  + e-^^at  + • > . 3.9 

By  linearity  of  the  DFT,  the  corresponding  sampled 
spectra  of  a sampled  segment  is  given  by  the.  sum  of 
two  terms  at  +f_ 

“ a 

eJa  -jJ(k-m-d) 

W(k)  - ijj-  e sindN(k-m-d) 

/ e-J«  lij(k+;r.+d> 

j + ■ *' j|'*"  e sind^(k+m+d)  . 3»10 

To  avoid  unnecessary  complication  in  the  math,  the 
following  arguments  will  be  in  terms  of  generating  a 
single  complex  exponential.  The  corresponding  real 
sinusoid  is  easily  obtained  by  summing  a pair  of 
conjugately  symmetric  exponentials  at  +f_,  the 

Cl 

desired  frequency. 

The  required  synthesis  algorithm  is  thus  to  add 
contributions  to  every  spectral  cell  (for  each  discrete 
frequency  desired)  according  to  3.8  (or  3.10).  Since 
there  will  be  exactly  m+d  cycles  of  the  required 
sinusoid  in  each  N sample  segment,  the  phase  angle 
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a must  be  Increased  by  2nd  radians  for  the  generation 
of  successive  segments  to  maintain  continuity  from 
one  segment  to  the  next. 

B.  An  M Cell  Approximation 

The  computation  of  the  W(k)  for  N/2  spectral  cells 
to  generate  N output  samples  obviously  requires  more 
computational  effort  than  direct  generation  of  the 
sinusoid  in  the  time  domain,  even  before  considering 
the  overhead  of  computing  the  inverse  DFT.  Thus  the 
exact  solution  using  FFT  techniques  is  not  practical 
unless  the  required  frequencies  can  be  constrained 
to  the  set  f&  ® nfQ,  n-0,  +1,  +N/2-1. 

If,  on  the  other  hand,  something  less  than  perfec- 
tion is  acceptable,  an  approximate  solution  may  yield 
the  hoped  for  computational  savings.  'Since  most  of 
the  power  in  the  spectral  domain  is  contained  in  cells 

near  the  desired  frequency,  it  would  seem  reasonable 

* 

to  approximate  the  complete  spectral  description  by 
truncating  the  tails  of  the  sind^Cx)  envelope  or  by 
some  other  equivalent  operation.  This  could  result  in 
an  acceptable  number  of  coefficients  to  be  specified 
for  each  segment  of  each  discrete  required  in  the 
output.  If  a 1024  complex  to  2048  real  IFFT  is  being 
utilized  and  if  only  ten  non-zero  coefficients  per 
discrete  will  produce  acceptable  output  results,  the 
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generation  of  10  complex  coefficients  in  the  spectral 
domain  becomes  quite  attractive  relative  to  having  to 
generate  2048  numbers  in.  the  time  domain  for  each 
sinusoid  desired. 

Assume  that  only  the  M largest  (in  magnitude) 
contributions  per1  complex  exponential  are  to  be 
retained  unmodified  in  the  sampled  frequency  domain. 
Since  the  contributions  of  all  other  cells  are  to  be 
eliminated,  the  total  power  is  decreased  by -the  sum 
of  the  squared  magnitudes  of  the  eliminated  spectral 
components.  In  fact,  due  to  the  orthogonality  of 
the  Fourier  frequencies,  the  M largest  components 
form  a least  mean  squared  error  estimate  of  the  original 
signal  if  only  M cells  may  be  non-zero.  Thus  if  mean 
squared  error  is  an  acceptable  criteria,  adding  the  M 
largest  contributions  based  on  either '3.8  (for  single 
sided  spectra)  or  3.10  '(for  double  sided  spectra) 
yields  the  optimum  approximation. 

Since  the  sindN(x)  envelope  of  the  spectral 

magnitudes  decreases  monotonically  on  either  side 

of  the  desired  frequency,  zeroing  all  but  the  M 

largest  components  is  equivalent  to  multiplying  the 

spectrum  by  a pulse  or  "window”  of  unit  magnitude 
' fs 

and  width  M-yj-  . Multiplication  in  the  frequency 
domain  is  equivalent  to  convolving  the  time  function 
(sinusoidal  segment)  with  the  transform  of  the  window. 
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The  resulting  time  function  segment  can  be  determined 
by  an  inverse  DFT  of  the  windowed  spectrum 


s(n)  * l W(k)  e 


,2irhk 


3.11 


k*m-^+l 


Using  3.6, 


s(n) 


-4 


Y2*  h 


I?  eJa  "y1  ' 1 f-T 


, ' M ,,  Z-0 

k«m-^+l 


3.12 


Interchanging  the  order  of  summation  and  refactoring  the 
exponentials, 


s(n)  * e 


i*  K;1 1 Jrt,t4)  Jt*"-'! 

Jo  * e J m ■,  e 


k=m-j+l 


3.13 


Factoring  an  e from  the  last  summation  and 

modifying  the  summation  limits  correspondingly, 


s(n)  « e 


£=0 


M 


N-l  1 J^(m+d)  j‘2g2L(n-t)  ? j2*jL(n-l) 

Z.  N e e Z„  e 


k— 1+1 


3. 1H 
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Factoring  the  original  sinusoid  from  the  summations,  and 
reordering  again. 


i«  j¥(m+d)  f J£r-(k_dli  N;1 

.A,  *ti  0 


• M 


Y2*n, 


Y2w , 


s(n) 


-d) 


3.15 


The  last  sum  is  Just  the  sind^(k-d),  so 


s(n) 


; m * 

. j«  eJ¥(m+d)  f j'¥(k-d>  -jjcH-iKk-d) 

* ew  e 2 e e 


k-^fl 


sindN(k-d) . 


3.16 


or,  finally,  the  desired  form, 


s(n)  = e^° 


j2jn(m+d)  I 

k=-|fl 


jjkk-d)(2n-N+l) 
e ' sindN(k-d). 


3.17 


This  is  just  the  desired  sinusoidal  segment,  but  modulated 
by  a gain  factor  which  is  dependent  on  both  n and  d. 

An  example  is  helpful  in  illustrating  the  result 
of  the  uniform  window  approximation.  Figure  3.1  shows 
a 128  sample  function  segment  which  covers  10.25  cycles 
of  a sinusoid.  Figure  3.2  is  the  corresponding  64  point 
single  sided  spectra  derived  by  DFT.  Note  the  180° 
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phase  difference  for  components  Just  above  and  Just 
below  the  actual  frequency.  A uniform  window  of  width 
M»10  is  applied  to  produce  the'  approximation  spectra 
in  Figure  3.3*  An  inverse  DFT  then  gives  the  approximation 
time  segment  shown  in  Figure  3.4. 

Crudely  speaking,  the  major  effect  of  truncating 
the  tails  of  the  segment  spectrum  is  to  attempt  to 
"match  up"  the  two  ends  of  the  segment.  This  is  to 
be  expected  since  the  "end-around"  discontinuity,  or 
equivalently,  the  discontinuity  from  one  "cycle"  to 
the  next  of  the  periodic  waveform  represented,  cannot 
be  sustained  without  the  frequency  components  contained 
in  the  tails  of  the  spectrum.  The  error  for  the  example 
given  is  shown  in  Figure  3*5,  and  is  representative 
of  the  relative  size  of  the  error  at  the  ends  of  a 
segment  versus  the  center  portion. 

A composite  time  segment  was  constructed  by  concatenating 
four  segment  approximations  derived  as  in  the  example 
above,  but  with  initial  phases  of  0,  ir/2,  it,  and  3ir/2 
radians,  respectively.  The  resultant  time  block  contains 
512  samples  covering  a total  of  4l  cycles  of  the  desired 
sinusoid,  and  is  shown  in  Figure  3.6.  Note  the  rather 
severe  effective  modulation  of  the  signal  envelope  in  the 
vicinities  of  the  segment  boundaries.  The  error  function 
for  the  composite  waveform  is  shown  in  Figure  3.7,  and 
has  a maximum  value  of  .55  relative  to  a perfect  (all 
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spectral  components  preserved)  waveform  amplitude  of  +1. 
The  maximum  excursion  of  the  approximation  waveform  has 
a relative  amplitude  of  1.09. 

The  spectrum  of  the  composite  approximation  wave- 
form is  shown  in  Figure  3.8.  The  largest  undesired 
component  is  about  2k  db  below  the  desired  signal 
strength.  Note  that  the  major  error  components  are 
at  frequencies  fa  + n/NT_,  where  f_  is  the  desired 

frequency  and  NT  is  the  length  of  each  approximation 

s 

segment.  This  is  the  behavior  expected  in  the  frequency 
domain  corresponding  to  the  effect  described  above  of 
the  time  waveform  being  effectively  modulated  on  a 
once  per  segment  rate. 

If  the  criterion  of  acceptability  is  other  than 

minimum  mean  square  error,  such  as  ratio  of  desired 

signal  to  largest  undesired  component  or  minimum 

amplitude  modulation  of  the  time  function,  the  uniform 

window  may  not  be  the  best  solution.  Shading  of  the 

♦ 

window  coefficients  was  briefly  pursued,  but  in 
general  decreased  side  lobe  strengths  were  obtainable 
only  at  the  expense  of  increases  in  amplitude 
modulation  in  the  time  function.  In  addition,  the 
fact  that  the  spectral  truncation  operation  reduces 
the  "end-around”  discontinuity  for  each  segment  means 
that  discontinuities  are  created  in  the  composite  output 
waveform  at  the  segment  boundaries.  Although  increased 
fidelity  can  be  achieved  by  using  a wider  window,  any 
computational  advantages  of  the  technique  are  soon  lost. 

k5 
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It  appears,  then,  that  this  particular  FFT  based  algorithm 
is  not  very  useful  unless  a rather  coarse  approximation 
is  acceptable,  preferably  with  only  a mean  square  error 
fidelity  requirement. 

C.  The  "Overlapped-Hanned"  Algorithm 

The  fact  that  the  major  difficulty  encountered  in 
the  above  algorithm  was  an  effective  modulation  of  each 
segment  suggests  the  possibility  of  purposely  introducing 
a specific  modulation  characteristic  and  to  compensate 
by  overlapping  and  summing  adjacent  segments;  i.e., 
introduce  redundant  processing.  This  is  all  the  more 
attractive  if  the  desired  discrete  components  are  to 
be  combined  with  a specified  broadband  signal,  since  it 
was  shown  in  Chapter  2 that  a two-times  redundancy  factor 
yielded  a fairly  satisfactory  broadband  algorithm.  If 
the  same  factor  of  2 can  be  exploited  for  the  discrete 
case,  the  approximation  coefficients  can  be  added  to 
the  broadband  spectra  before  transforming,  and  thus  one 
IFFT  does  all. 

Since  the  contributions  in  each  of  the  overlapped 
segments  add  coherently  to  form  the  total  signal 
for  the  discrete  case,  the  sine  pulse  weighting  used 
for  the  broadband  algorithm  is  not  applicable.  A 
triangular  pulse  modulation  form  meets  the  requirement 
in  the  time  domain,  but  is  not  obtainable  by  simple 


47 


NOLTR  74-215 


manipulation  of  the  spectrum  of  the  signal.  On 
the  other  hand,  the  raised  cosine  modulation  (corresponding 
to  "Hanning"  shading)  of  the  segment  appears  to  be  an 
almost  perfect  solution  to  the  problem.  Overlapped 
signal  segments  that  have  been  modulated  by  the  Hanning 
pulse  pH(t) 


.5  + .5  cos  2irt/NT  ; |t|  < NT/2 
o;  i 1 1 1 nts/2 


3.18 


will  add  perfectly  to  produce  the  original  waveform. 

The  question  thus  becomes:  can  the  Hanned  segment  be 

more  easily  approximated  than  the  unmodulated  segments? 

Returning  to  the  spectrum  of  the  unweighted  segment. 
Equation  3.8,  the  spectrum  of  the  Hanned  segment  WH(k) 
is  obtained  by  convolving  with  the  transform  of 
^•(1  + cos  2ir  n/N).  The  latter  is  the  well  known 

pH<>0, 


PH(k)  » jj-  S(k+1)  + | «(k)  + «(k-l)  3.19 

i.e.,  impulses  at  zero  frequency  and  at  + one 

fs 

frequency  cell  -vr  • 

Convolving  3.8  with  3.19  the  result  is 
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1 


\ 

i 

I 

j 


WH(k)  * eJ  e 


j«.  J,(-1 TT)d 


■ t (k-m) 

(i*  e w sind^ (k-m-d) 


l “Jir(^)(k-.m+l) 


+ te 


sindjj(k-l-m-d) 


sindjj(k+l-m-d)+^-  e 


1 -j*  (k-m-1) 


3.20 


Manipulating  and  writing  out  the  sind  functions, 

,N-1^  _,_/N-l 

i < „ J w ' 

_ f’lr\  = ® P 

H 


1 la  J (“TT”^  “jir  (-rj— ) (k-m) 
Wu(k)  = ^ eJa  e N e N 


slnir  (k-m-d) 


s in  jj-(  k-m-d) 


. 1 sin[.Ck-m-d)+ir]  , 1 .^TT5 

+ e r — — + *•  e 

''  sin[jj(k-m-d+l)] 


rN— 1 < 


slnCff  (k-m-d)-Tr] 


sinC^(k-m-d-l)] 


3.21 


In  the  last  two  terms,  sin[ir (k-m-d )+ir]  is  just 
-sinO  (k-m-d)],  which  can  then  be  factored  out  of 
the  brackets.  Since  the  approximation  is  concerned 
with  the  WR(k)  for  k near  m,  sin  tr— is 
approximately  equal  to  n for  reasonably  large 

N,  say  on  the  order  of  10  . Similarly  the  factor 

appearing  in  the  exponentials  is  approximately 
unity  for  N large.  Thus,  for  k in  the  vicinity  of  m, 
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WH0C)  : £ eJ‘ 

+ 2 (k-m-d+1 ) + 2tk-ra-a-l) } 3.22 

Factoring  a from  inside  the  brackets  and  combining, 

w fki  * l eJ  (•+»<*).  ( t n (k~m)/ sinu  (k-m-d)  , f2(k-m-d)2-l}  j 
wHu;  ^ ? e k i)  \ „ ck-m-d)  v (k-m-d)2-l  > 

3.23 

The  magnitudes  for  the  frequency  contributions  around  k = m 
are  given  by  the  factor  in  brackets.  The  phases  are  determined 
by  a + itd,  plus  some  number  of  it  rotations  determined  by 
k - m in  the  last  exponential,  and  the  sign  of  the  factor 
in  brackets.  In  the  example  given  below,  it  is  shown 
that  these  factors  of  (-l)n  are  responsible  for  three 
consecutive  phase  shifts  of  it  radians  each  in  the  vicinity 
of  k = m. 

4 

The  approximate  WH(k)  given  by  3.23  to  be  generated 
for  the  synthesis  algorithm  has  the  following  exploitable 
characteristic.  As  was  noted  before,  only  the  phase  a 
changes  for  any  of  the  W^(k)  from  one  segment  to  the  next. 
Furthermore,  the  only  phases  required  are  either  o + itd  or 
a + itd  + it.  Thus,  if  M magnitudes  are  computed  and  stored 
for  an  M cell  approximation,  only  one  cosine  and  one  sine 
are  needed  to  determine  the  real  and  imaginary  parts  to 
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be  added  to  the  synthesised  spectrum  for  each  segment. 

Note  that  unless  a particular  phase  is  required,  the  initial 
phase  (a+ird)  for  the  first  segment  can  be  arbitrary,  and 
all  that  is  required  thereafter  is  to  increase  the  a+ird 
term  by  ird  on  each  iteration. 

Since  the  primary  virture  of  Hanning  weighting  is  to 
greatly  reduce  the  sidelobes  of  the  sinc(x)  DFT  filter 
function,  it  provides  almost  the  same  reduction  in  the  tails 
of  the  sind^Cx)  form  of  3.8.  The  result  is  that,  although 
as  before  perfect  synthesis  requires  contributions  in  all 
frequency  cells, the  Hanning  modulated  discrete  segment 
achieves  a very  close  approximation  with  many  fewer  non-zero 
components . 

The  basic  parameters  of  the  previous  example  were 
repeated  in  the  example  below  illustrating  the  application 
of  the  "overlap-Hanned"  algorithm.  A set  of  time  function 
segments,  each  consisting  of  10.25  cycles  in  128  points, 
with  phases  incremented  by  2tt/8  rad.  for  each  successive 
segment,  are  forward  transformed,  Hanned,  tails  truncated, 
inverse  transformed  and  combined  to  form  a composite  output 
time  function  of  4l  cycles  in  512  points  as  before.  Figure 
3.9  shows  the  same  10.25  cycle  segment  as  shown  in  3.1,  but 
after  modulation  by  the  raised  cosine  pulse,  of  "Hanning’1 , 

The  corresponding  spectrum  is  shown  in  Figure  3.10.  Note 
that  the  tails  of  the  spectrum  have  been  reduced  considerably, 
such  that  barely  four  components  on  each  side  of  the  actual 
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frequency  are  within  the  48  db  range  plotted.  Note  in 
addition  that  the  phase  undergoes  three  shifts  of  if 
radians  each  in  the  immediate  vicinity  of  the  desired 
frequency,  compared  to  the  single  it  radian  shift 
(Figure  3.2)  for  the  unHanned  segment. 

The  tails  of  the  spectrum  of  each  Manned  segment  are 

then  truncated  to  produce  the  spectrum  of  the  approximation 

as  before.  For  this  example,  only  the  six  largest  components 

in  Figure  3.9  are  retained,  which  corresponds  to  throwing 

away  all  contributions  that  are  more  than  about  40  db 

below  the  largest  one.  It  is  interesting  to  note  that, 

in  Figure  3.2,  no  components  in  the  unHanned  segment 

are  below  the  -40  db  level.  Finally  Figure  3.11  shows 

the  composite  512  sample  approximation  waveform  obtained  by 

the  overlap-Hanned  algorithm.  The  error  is  plotted  in 

Figure  3.12.  The  approximation  spectrum,  shown  in 

Figure  3.13  has  no  sidelobe  components  larger  than  50  db 

below  the  desired  component.  Figure  3.14  is  the  spectrum 

of  the  error  alone,  which  consists  primarily  of  components 

displaced  by  multiples  of  2/NT  , or  the  reciprocal  of  half 

s 

the  segment  length.  The  doubling  of  the  effective  modula- 
tion rate  relative  to  the  first  algorithm  is  due  to  the 
factor  of  2 overlap,  since  the  modulation  of  each  segment 
is  definitely  still  characterized  by  the  period  NT  . 

The  "overlap-Hanned”  synthesis  algorithm  for  discrete 
components  yields  very  high  fidelity  with  a very  small- 
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number  of  non-zero  spectral  contributions 
a six  component  approximation  empirically 
a maximum  error  of  slightly  more  than  If, 
spectral  sidelobes  of  about  -50  db.  Even 
approximation  has  less  than  5f  error,  and 
max  sidelobe  levels. 


0 


In  particular, 
appears  to  yield 
and  maximum 
a 4 cell 
about  -35  db 
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CHAPTER  IV 

TIME  DOMAIN  SYNTHESIS  OF  DISCRETE  FREQUENCIES 
A.  Harmonic  Sets 


Quite  often  it  is  desirable  to  generate  sets  of 
narrowband  components  which  are  harmonically  related.  That 
is,  all  of  the  members  of  a particular  set  are  Just  the 
spectral  components  Yk  of  a periodic  or  almost  periodic 
time  waveform.  Let  the  tine  function  of  a particular 
harmonic  set  be  denoted  by  y(t),  defined  by 


y(t) 


k»-*» 


J2ikf0t 


4.1 


where  fQ  is  the  fundamental  frequency  of  y(t).  The 
spectrum  of  y(t)  is  then 


Y(f)  * l Yk  6 (f  - kf0)  . 4.2 

k=-» 

Although  each  of  the  components  of  Y(f)  could  be 
generated  explicitly  by  the  techniques  described  in 
Chapter  3,  it  becomes  profitable  to  generate  y(t) 
directly  in  the  time  domain  if  the  number  of  non-zero  spectral 
components  Y(kfQ)  is  fairly  large. 
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If  identical  copies  of  a waveform  defined  over  an 
interval  T are  concatenated  in  time,  the  result  is  a 
periodic  signal  with  fundamental  frequency  1/TQ. 

Thus  if  a lossless  delay  line  of  length  TQ  containing 
one  period  of  a desired  waveform  were  recirculated  as 
illustrated  in  Figure  4.1  the  signal  observable  at  any 
tap  would  be  the  periodic  signal  desired,  with  fundamental 
frequency  f * 1/TQ.  If  the  waveform  contained  in  the 
delay  line  is  exactly  one  full  cycle  of  a sinusoid,  the 
spectrum  of  the  signal  generated  will  be  non-zero  only 
at  the  fundamental  frequency  +f  . However,  if  an 
arbitrary  waveform  is  contained  in  the  delay  line,  then 
in  genera.!  the  spectrum  may  contain  any  of  the 
frequencies  harmonically  related  to  fQ,  i.e»,  kfQ  for 
all  integer  k.  The  spectral  values  at  these  harmonic 
frequencies  are  just  the  coefficients  obtained  by 
expanding  the  TQ  seconds  of  signal  in  the  delay  line 
in  a Fourier  series. 

The  discrete  equivalent  cf  the  above  operation  is 
obtained  by  loading  a shift  register  with  the  samples  of 
one  period  of  a waveform  y(t),  connecting  the  last  stage 
to  the  input  for  recirculation,  and  observing  a discrete 
periodic  signal  yx(t)  at  any  stage  in  the  shift  register. 
For  an  N stage  register  and  a shift  rate  fr,  the  sampled 
periodic  signal  produced  has  a fundamental  frequency  f » 
f^/N.  Since  the  resultant  waveform  is  both  sampled  and 
periodic,  its  spectrum  Yx(f)  is  also  both  sampled  and 
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CONTENTS  OF  LINE  AT  TIME  t) 

OR  ONE  PERIOD  OF  OUTPUT  SIGNAL 


FIG.  4.1  GENERATION  OF  PERIODIC  SIGNAL  BY  RECIRCULATING  DELAY  LINE 
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periodic,  with  a sample  spacing  of  f and  period 
NfQ  = fr>  Since  a shift  register  may  be  shifted 
at  a variety  of  rates,  the  fundamental  frequency  of  the 
signal  generated  by  the  discrete  system  may  be  easily 
adjusted  to  any  desired  value  by  the  appropriate 
choice  of  shift  rate  fp.  Figure  *1.2  illustrates  the 
output  of  a six  state  register  for  two  different 
shift  rates,  and  the  corresponding  spectra  for  each. 

If  the  shift  register  version  described  above  is 
applied  to  the  generation  of  a periodic  signet  with  an 
arbitrary  harmonic  structure,  care  must  be  taken  that  the 
shift  register  is  sufficiently  long  to  allow  independent 
definition  of  all  desired  harmonics.  Since  the  discrete 
spectrum  Yx(f)  is  periodic  with  period  NfQ,  it  is 
obvious  that  there  are  at  most  N different  complex 
harmonic  amplitudes.  However,  since  yx(t)  is  a real 
time  series,  the  spectral  components  exhibit  conjugate 
symmetry  about  f = 0;  i.e.,  if 

j2irkt 

„ *-ijr — 

yx(t)  = l Yk  e 0 = Real  4.3 

k~-» 


then 


4.4 
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(2b) 


FIG.  4.2  OUTPUT  OF  6 STAGE  SHIFT  REGISTER  FOR  TWO  DIFFERENT 
SHIFT  RATES 
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where  * denotes  the  complex  conjugate.  Thus  for  a N 
stage  shift  register,  at  most  [(N-l)/2]  harmonics  plus 
the  d.c.  component  are  independently  specifiable,  where  [ ] 
denotes  the  integer  part. 

A further  complication  exists  when  the  valid  output 
band  may  be  significantly  wider  than  the  band  coye'ed  by 
the  desired  harmonics  of  the  generated  signal. 

Figure  4.3  illustrates  a ease  where,  if  zero  filled 
harmonics  are  not  explicitly  included  in  the  waveform  stored 
in  the  shift  register  (corresponding  to  sampling  at  greater 
than  the  Nyquist  rate  and  therefore  additional  len-^h  of  the 
shift  register),  undesired  harmonics  are  created  ii  the 
output  band  of  interest  for  low  shift  rates  (low  fundamental 
frequencies)  due  to  the  periodic  nature  of  the  generated 
signal  spectrum.  In  general,  if  the  required  valid  output 
bandwidth  is  B and  the  lowest  fundamental  frequency  required 
is  f the  minimum  number  of  harmonics  that  must  be 
specified  is  B/f^in,  or  equivalently,  the  minimum  shift 
register  length  is  2B/fmin»  From  Figure  4.3c  it  can  be 
seen  that  if  the  highest  non-zero  desired  harmonic  is  m, 
then  the  minimum  value  the  shift  rate  can  be  is  B + mfmin, 
or  the  valid  output  band  plus  room  for  the  image  spectra  at 
the  shift  rate  f^.  The  minimum  length  of  the  shift 
register  is  then  given  by  the  next  integer  larger  than 

VW  J-o-  Nmin  ■ >VW+1  ■ tB/W**1- 

Thus,  for  all  fundamental  frequencies  for  which  the 
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desired  waveform  must  be  generated,  the  unwanted 
harmonics  are  outside  the  band  of  interest  and  may  be 
removed  by  analog  filtering  after  D/A  conversion,  or 
just  ignored. 

B.  The  "Slip-Sample”  Algorithm 

The  above  technique  for  generating  periodic  signals 
with  arbitrary  harmonic  structure  and  variable  fundamental 
frequency  requires  a dedicated  recirculating  shift 
register  and  a controllable  rate  shift  clock-,  both  of 
which  are  special  purpose  hardware.  The  following 
technique  is  an  adaptation  of  the  principle  of  the 
above  implementation,  but  is  designed  to  be  realized 
as  software  in  a general  purpose  computer. 

Assume  that  the  required  number  of  discrete  samples  of 
one  period  of  the  desired  signal  is  stored  in  a block  o< 
random  access  computer  memory.  If  the  samples  are 
sequentially  read  from  the  block  by  the  computer  and 
fed  to  the  output  D/A  converters,  the  result  is  identical 
to  the  shift  register  output  described  above  (the  computer 
returns  to  the  beginning  of  the  block  after  each  pass  through 
the  block).  The  rate  at  which  the  samples  are  read  and 
output  is  determined  by  a clock  interrupt  or  other 
timing  technique,  and  as  before  determines  the  fundamental 
frequency  of  the  output  signal. 
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In  a more  general  case,  however,  the  generation  of 
the  periodic  signal  may  be  only  a part  of  the  overall 
signal  processing  problem,  and  may  be  fed  directly 
to  other  parts  of  the  sampled  data  system.  In  the 
usual  case  the  entire  system  will  be  operating  at 
one  fixed  sample  rate  f_,  and  each  part  of  the  system 
must  accept  and  process  data  at  the  rate  fg.  This 
rate  is  not  in  general  related  to  the  desired  fundamental 
frequency  f of  the  periodic  signal  to  be  generated. 

One  possible  solution  to  this  constrained  problem 
would  be  to  adjust  the  size  of  the  memory  block  to 
be  such  that  at  the  sample  rate  f , once  through  the 
block  would  be  TQ  = l/fQ  seconds,  where  f is  the 
desired  fundamental.  This  has  two  rather  severe 
drawbacks.  First,  the  realizable  fundamental  frequencies 
are  given  by  f /n,  where  n is  an  integer.  Thus,  except 
for  extremely  low  fundamental  frequencies  (large  n), 
the  realizable  frequencies  are  few  and  far  between. 

In  addition,  even  for  very  low  fundamentals  where 
reasonable  frequency  resolution  is  available-,  each 
< e in  fundamental  frequency  requires  rederiving 
the  complete  set  of  time  samples  stored  in  the 
memory  block. 

A more  viable  solution  to  the  fixed  sample  rate 
problem  is  to  maintain  a fixed  block  size  adequate 
for  the  lowest  desired  fundamental  frequency,  and 
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to  vary  the  effective  rate  at  which  the  computer 
scans  through  the  block.  This  variation  in  readout 
rate  is  accomplished  by  maintaining  and  incrementing 
the  position  in  the  memory  block  with  a precision 
much  greater  than  one  memory  cell.  Since  it  is  not 
really  possible  to  read  between  memory  cells  or 
desirable  for  computational  loading  considerations 
to  interpolate,  the  block  position  is  rounded  (or 
truncated)  to  the  nearest  integer  memory  location  for 
the  actual  data  access.  In  simple  terms,  the  net 
effect  is  to  "stretch"  the  stored,  waveform  by 
occasionally  duplicating  a value  in  the  output  (reading 
a memory  cell  twice  before  going  on)  or  to  "shrink" 
the  waveform  by  occasionally  skipping  a cell  in  the 
readout . 

Consider  a memory  block  of  length  N and  block 
position  parameter  P,  as  shown  in  Figure  4.4.  Assume 
that  P is  maintained  to  a sufficiently  high  resolution 
as  to  be  essentially  continuous  relative  to  the 
quantization  of  the  memory  into  N words.  If  the  block 
of  memory  contains  one  period  of  some  waveform,  then 
P is  in  effect  a phase  angle  of  the  fundamental  of  the 
periodic  waveform.  For  a non-integer  P between  the 
4th  and  5th  samples  as  Illustrated,  the  contents 
of  the  fourth  cell  are  used  as  the  output  for,  say, 
the  ivn  sample.  P is  then  incremented  by  an  amount 
F derived  from  the  desired  fundamental  frequency. 
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MEMORY  CELL  CORRESPONDING  TO  SAMPLE  AT 
TIME  4 ACCESSED  FOR  4 £ P < 5 


0 1 2 3 4 5 N-2  N-1  N 


i t t t m i i i ;i 

P = 4.437  . . . (IN  UNITS  OF  MEMORY  LOCATIONS) 


( <(>  » 2n  . 4.437  ...  (AS  A PHASE  ANGLE  IN  RADIANS) 

N * 

I 

l 


V 


J 


VALID  RANGE  OF  P 


FIG.  4.4  SCHEMATIC  REPRESENTATION  OF  CONTINUOUS  POSITION  IN  MEMORY 
BLOCK  AS  PHASE  ANGLE 
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and  the  appropriate  memory  cell  is  accessed  for  the 
(i+1)  sample.  P is  specified  with  the  same  resolution 
as  P,  and  is  therefore  also  not  in  general  an  integral 
number  of  memory  cells.  After  incrementing  by  F,  P 
is  retained  modulo  N to  maintain  it  within  its  valid 
range.  The  fundamental  frequency  f is  then  given 
by  the  rate  of  passing  completely  through  the  b3ock, 
or 


It  is  important  to  note  here  that  the  memory 
locations  accessed  twice  (or  skipped,  as  the  case  may 
be)  are  in  general  different  on  each  pass  through  the 
block.  In  addition  the  number  of  "doubles”  or  "skips” 
per  pass  may  in  general  alternate  between  two  values 
such  that  the  average  number  per  pass  is  non-integral. 
Thus  this  technique  differs  significantly  from  just 
modifying  the  block  by  duplicating  or  dropping  one 
or  more  samples.  It  is  also  obvious  that  for  frequencies 
which  require  a non-integral  number  of  "glitches"  per 
pass,  the  output  signal  is  not  truly  periodic,  nor  is 
it  even  composed  of  samples  of  the  desired  periodic 
signal.  It  is  really  only  an  approximation  to  the 
sampled  version  of  the  desired  signal  and  on  the  average 
(in  some  sense)  has  the  right  behavior. 
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The  analysis  of  the  sampled  data  signal  generated 

by  this  technique  is  aided  by  identifying  the  operations 

performed  as  components  of  a more  familiar  system. 

Let  the  desired  periodic  signal  be  y(t),  with  period 

TQ  and  fundamental  frequency  f = 1/T0.  Assume  that 

samples  are  taken  from  y(t)  at  a rate  f which  is 

S1 

N times  the  fundamental.  Thus  there  will  be  exactly 

N samples  in  each  period  of  the  sampled  signal  yx(t). 

The  samples  (or  impulses)  of  y (t)  are  now  applied 

to  the  input  of  a "boxcar-integrator"  filter,  i.e., 

a filter  with  an  impulse  response  of  a rectangular 

pulse  of  unit  amplitude  and  duration  TQ  centered  at 

t = 0.  Let  the  filter  output  be  designated  yfa(t). 

If  yb(t)  is  now  resampled  at  a second  rate  fg  , in 

general  not  related  to  f , samples  taken  during  any 

s - 

« r 

interval  (n-^)TQ  < t < (n+^OT  will  yield  the  value 
of  the  sample  taken  'from  y(t)  at  t = nTQ.  This  is 
exactly  analogous, to  the  situation  where  samples 
taken  from  y(t)  are  stored  in  a computer  memory,  and 
then  accessed  by  rounding  sample  "times"  to  the 
nearest  memory  address,  assuming  here  that  samples 
are  stored  such  that  increasing  time  corresponds  to 
Increasing  addresses,  (Truncating  instead  of 
rounding  to  form  the  address  would  introduce  a delay 
of  half  a memory  cell,  but  would  yield  equivalent 
overall  results.)  The  output  signal  is  therefore 
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the  result  of  1)  multiplying  the  desired  y(t)  by  an 

impulse  train  of  frequency  f , 2)  convolving  the 

S1 

result  with  the  impulse  response  of  the  boxcar  filter, 

and  3)  resampling  its  output  by  multiplying  by  the 

impulse  train  of  frequency  f . Correspondingly, 

s2 

the  spectrum  of  the  output  signal  is  obtained  from  the 

spectrum  of  y(t)  by  1)  convolving  with  an  impulse 

train  of  spacing  f (which  is  some  integral  multiple 

51 

of  the  fundamental  frequency)  as  a result  of  the  initial 

sampling  process,  2)  multiplying  by  sinc(f/f  ), 

S1  1 

the  filter  function  of  the  boxcar  integrator,  and 

3)  convolving  the  result  with  a second  impulse  train 

of  spacing  f . 

s2 

Since  for  the  case  of  interest  here  the  signal  y(t) 
is  periodic,  its  spectrum  Y(f)  is  impulsive  with 
spacing  fQ: 


Y(f) 


N/2-1 

' I 

k=-N/2+l 


Yk  5(f  - kf0> 


4.6 


where  Y(f)  has  been  explicitly  band-limited  in 
anticipation  of  sampling  at  Nf  . Performing  the 
sampling  operation,  Y(f)  must  be  duplicated  every 
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Yx(f)  * l Ytt-raf  ) 
m»—  S1 


- N/2-1 

l I Y.  Ck+mN])  4.7 

m«—  k=-N/2+l  * 0 


Multiplying  by  the  filter  function. 


, * N/2-1 

Y (f)  = 1 SincCf/Nf0)  I I Vk  «(f-f0[k+mN])  4.8 

o m=-»  k=-n/2t1 


Convolving  with  the  second  impulse  train  for  the  final 
step,  the  spectrum  of  the  output  y£(t)  is  obtained: 

Yj(f)  = l Y.  (f-rf  ) 4.9 

D r=—  D s2 


or 


• ’ f-rfs 

1 sincO™-.2--)  I 

o o m=- 


N/2-1 

l , 

k— N/2+1 


Yk  6 ( [f-rf,  ] 

K 02 


- [k+mN]f0)  . • 4.10 

The  system  under  consideration  has  a fixed  system  sample 

rate  f , and  thus  the  second  sample  rate  f is  constrained 
s s2 

to  the  value  f = f . Now  define  a parameter  a as  the 

s 

ratio  of  the  first  sampling  rate  f to  f 
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nr 


. af  s 

“TT 


4.11 


Substituting  4.11  into  4.10, 


„ « , f-rfe  « N/2-1 

Y>(f)  ‘ rL  ^ slncC"^f)  mL  J.N/2+1 


- Cg+m]afs) 


4.12 


The  resultant  spectrum  Y^(f)  has  a replica  of  the  boxcarred 
spectrum  at  each  multiple  of  fg,  as  is  expected  for  any 
sampling  operation,  but  with  the  effective  frequency 
scale  of  each  replica  about  the  appropriate  value  of  rf 

s 

variable  with  a.  Since  a is  proportional  to  fQ,  the 

frequency  of  the  fundamental  of  the  harmonic  set  generated 

can  be  controlled  by  controlling  a.  From  4.11  f = of  , 

S<i  s 

IX  1 

or  jr-  - ci jr — , so  that  the  spacing  between  samples  for 

S S j / 

the  second  Sampling  process  is  a times  the  spacing  of 

samples  going  into  the  boxcar  filter.  The  latter  samples 

are  stored  one  to  each  memory  cell  in  the  computer  memory, 

and  thus  1/f  in  terms  of  memory  cell  spacing  is  unity. 

S1 

The  spacing  in  memory  cells  tetween  samples  for  the 
output  sampling  process  is  therefore  o.  But  this  spacing 
is  just  what  was  previously  defined  as  the  quantity  F, 
and  from  4.5  and  4.11,  F = a is  easily  verified. 
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The  initial  sampling  process  results  in  a non-band 

limited  spectrum  even  though  the  original  spectrum  Y(f) 

was  band-limited.  The  effective  boxcar  filtering 

operation  reduces  the  amplitude  of  the  replicas  generated 

by  the  first  sampling  operation,  but  since  the  sinc(f/f  ) 

S1 

envelope  only  falls  off  as  1/f,  considerable  power  still 

remains  in  the  undesired  replicas.  After  the  second 

sampling  operation  these  reduced  but  non-zero  replicas 

of  the  desired  spectrum  in  the  vicinity  of  +f  , +2f  , 

“ 81  " S1 

etc.,  of  the  box-carred  spectrum  may  be  aliased  back 

into  the  band  of  interest  by  the  resampling  operation, 

and  in  general  will  appear  as  weak,  discrete  frequencies 

not  har&*v  .ically  related  to  the  desired  fundamental 

tQ.  These  artifacts  may  be  quite  large  if  generated 

by  harmonics  in  the  original  waveform  with  frequencies 

close  to  1/2  of  the  original  sample  rate  f . This  is 

S1 

due  to  the  fact  that  the  sinc(f/f  ) envelope  of  the  box- 

S1 

carred  spectrum  is  down  only  about  4db  at  the  Nyquist 

frequency  f /2.  If  the  original  waveform  is  oversampled 
S1 

by  a factor  of  two  such  that  all  non-zero  harmonics  are 

below  f /4,  the  minimum  attenuation  provided  by  the 
S1 

sine  envelope  becomes  a little  better  than  lOdb,  and 
increases  by  6db  for  each  additional  factor  of  2 in 
oversampling  (see  Figure  4.5).  Thus  by  increasing  the 
memory  required  to  store  the  waveform  it  is  possible 
to  make  artifacts  as  small  as  desired.  However,  the 

6db  per  octave  roll-off  of  the  sine  function  makes 
this  a fairly  expensive  solution  to  the  problem. 


r3 


NOLTR  74-215 


0 r 


-10 


-20 


-30 


-40 


• -3.92 


• -10.45 


• -17.1 


• -23.6 


• -29.8 


• -36 


16 


-42  • 


J 


32 


64 

OVER- 
SAMPLE 
RATIO  r 


FIG.  4.5  MINIMUM  ATTENUATION  OF  ALIASED  SPECTRA  FOR  VARIOUS  OVER-SAMPLE 

RATES.  (f$)  = r • k • fo  ; WHERE  k • fa  IS  THE  HIGHEST  NON  ZERO  HARMONIC  DESIRED) 


If  the  memory  address  increment  F which  determines 
the  fundamental  frequency  of  the  output  is  fluctuated 
rapidly  by  some  small  amount,  t\a  effect  is  to  smear  ' 
the  power  at  each  of  the  harmonically  related  disc,  'ete 
frequencies  over  some  finite  bandwidth.  The  bandwidth 
of  the  ntn  harmonic  will  be  n times  the  bandwidth  of 
the  fundamental.  In  addition  to  being  useful  to  generate 
finite  bandwidth  narrowband  processes,  this  phenomena 
also  tends  to  smear  out  the  discretes  in  the  band  of 
interest  generated  by  aliasing.  Unfortunately  if  highly 
stable  discretes  are  desired,  the  aliased  contributions 
may  still  be  unacceptably  large. 

Mora  exotic  interpolation  schemes  or  other  prefiltering 
operations  could  be  performed  on  the  stored  waveform 
prior  to  resampling  to  reduce  aliasing,  but  even  linear 
interpolation  requires  a multiplication  per  sample,  and 
therefore  it  as  well  as  more  sophisticated  filtering 
schemes  are  impractical  to  Implement. 

C,  Stochastic  Filtering 

Although  conventional  filtering  algorithms  are 
impractical,  an  effective  prefiltering  operation  can 
be  performed  by  "jittered"  resampling,  and  thus  the 
aliasing  problem  can  be  reduced.  The  "jittered"  sampling 
technique  is  quite  attractive  relative  to  more  straight- 
forward prefi] tering  methods  due  to  its  extremely  simple 
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implementation.  Practical  application  to  the  Memory- 
Read-Out  signal  generation  problem  will  be  discussed 
after  deriving  the  spectrum  of  the  output  if  ’’Jittered" 
sampling  is  employed. 

The  problem  is  as  follows:  the  function  s(t)is 

given  with  Power  Spectral  Density  (PSD)  $ (f)  and 

s 

autocorrelation  function  (ACF)  (t).  Therefore 

s 


*s(t)  = 


4.13 


= /“Mt)  e"j27rfTdr 


4.14 


is  just  the  Fourier  Transform  of  4>_(t). 

s 

Samples  are  now  taken  of  s(t)  according  to  the  ’’jittered’’ 

t*  H 

sampling  scheme;  i.e.,  the  i sample  is  taken  at  the  time 
i T + v . , where  T = 1/f  is  the  period  associated  with 
the  underlying  basic  sample  rate  fs,  and  the  v1  ar« 
independent,  identically  distributed  random  variables 
with  probability  density  function  (PDF) 


Pv  (t)  = p(t)  , all  i 


4.1? 


The  sample  values  of  s(t)  obtained  are  then  assigned 
to  a regularly  spaced  impulse  train  of  period  T„  to 

o 

y 

form  a sampled  function  s (t)  given  by 


< sJL- : 
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sX(t)  = y s(nT  + v ) 6 (t  - nT  ) 4,i _b 

- n=-«  s “ 

The  problem  thus  is  to  determine  the  PSD  of  s (t)  in 

terms  of  * (f),  the  PSD  of  s(t);  and  of  p(t),  the  PDF  of 
s 

the 

Consider  the  sequence  autocorrelation  of  the  original 
"jittered"  samples  of  s(t): 

R(r)  = s(n!l?s  + vn)s((n+r)i’s  + vn+r)  4.17 

where  the  average  is  to  be  taken  over  all  n and  over  the 
vi.  For  vn  and  vn+r  identically  equal  to  zero,  R(rj  is 
Just  the  autocorrelation  of  s(t)  for  rTg  delay,  i.e., 

R(r)  = *s(rTg);  all  v 5 0 4.18 

For  an  arbitrary  distribution  function  p(t)  for  the 

4 

v . , R(r)  is  the  statistical  average  of  <j>  (t)  taken  over 

all  possible  vn+r  - vn.  For  r^O  the  are  independent, 

and  the  PDF  of  the  difference  of  two  of  the  is  just 

the  auto-correlation  of  the  PDF  of  each.  Let  pd(t)  be 

the  PDF  of  v , - v . i.e., 

„n+r  >n’ 

co 

P^(t)  = / p(a)  p(a+t)da  4.19 
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Then  the  desired  sequence  autocorrelation  function  R(r) 
is  given  by: 

CO 

R(r ) = / pd(t)  ^s(rTs  + t )dt ; r^O  . 4.20 

Since  p^Ct)  is  real,  positive  and  symmetric,  4.20  may 
equivalently  bo  written 

00 

R(r)  = / pd(t)  4>s (rTg  - t)dt;  r^O  . 4.21 


I 


If  the  values  of  R(r)  are  now  assigned  to  a regularl" 

spaced  impulse  train  of  period  T , the  discrete 

s 

autocorrelation  function  of  the  sampled  function  sx(t) 
is  obtained.  Let  Rx(t)  be  the  autocorrelation  function 
of  sx(t)  . Then 


RX(t)  = l 


r=*-» 

r^O 


R(r) 


6 Ct  - rTg ) + Rx (0 ) . 


4.22 


Rx(0)  is  unaffected  by  the  jittering  operation  and 
therefore  must  be  handled  as  a separate  case.  Substituting 
for  R(r) 


Pd(tHs(rT3-t)dt]5(T-rTs)  + Rx(0) 

4.23 
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Now  define  a continuous  autocorrelation  function 
R(t)  as  follows: 


R(t)  = /"  PdCt)4  (r-t)dt  . 4.24 

—m  U S 

If  R(t)  is  sampled  at  points  rT_,  the  result  is  Rx(t) 

s 

except  for  r=o.  Let  the  Fourier  transform  of  R(t) 
be  *R(f),  and  the  Fourier  transform  of  Pd(t)  be  P^Cf). 

Then  from  4.24 


*R(f)  = Pd(f)»s(f)  . 4.25 

Finally,  since  Rx(x)  is  (except  for  t=o)  just  a sampled 
version  of  R(t),  the  Fourier  transform  of  Rx(t)  is  given 
by 


v(f> 


Z + K * K + Z ZtU-rtJ*  U-rr  ) 

PC««P9 


4.26 


The  constant  K arises  due  to  the  fact  that  Rx(0)  / R(0) 
since 


Rx(0)  = E{[S(nTs  + v±)l2}  « *gC0)  4.27 
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and 


R(0)  = I"  PflCt )♦  (t )dt  4.28 

— Oft 

which  is  equal  to  $g(0)  only  if  pd(t)  ^(t).  Thus  R (t) 
can  be  written 

Rx(t ) = l R(t)  6(T-rTs)  + [*s CO ) - RCO)]  6(x).  4.29 

J»=  — CO 

The  transform  of  the  last  term  in  4.29  yields  the  constant 
K 


K - *s(0)  - R(0)  , 4.30 

or,  using  4.28, 

00 

K « $s(0)  - / pdCt)  *s(t)  dt  . 4..31 

— Oft 

♦ 

The  effect  of  the  "jittered”  sampling  can  be  seen 
from  4.25  and  4.26  to  be  almost  equivalent  to  passing 
s(t)  through  a filter  with  an  impulse  response  equal 
to  the  PDF  of  the  v^,  p(t).  However,  the  one  anomaly 
of  stochastic  filtering  is  the  white  noise  term  K 
which  is  added  to  the  PSD  of  the  final  output.  The 
significance  of  the  white  noise  term  is  apparent  if 
one  realizes  that  the  stochastic  filtering  is  a power 
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lossless  operation,  and  that  the  power  missing  from  the 
attenuated  frequencies  shows  up  as  white  noise.  Since 
the  objective  was  to  reduce  the  discretes  generated  in 
the  band  of  interest  by  aliasing,  trading  them  off  for 
white  noise  is  usually  quite  acceptable.  Further 
properties  of  '’jittered1'  sampling  are  discussed  in 
Balakirshnan^  and  Shapiro  and  Silverman^. 

The  implementation  of  the  stochastic  filter  by 
jittered  sampling  in  the  memory  read-out  algorithm 
is  fairly  simple  compared  with  the  computational 
requirement  of  normal  filtering  methods,  assuming  that 
the  are  readily  available.  After  the  block  position 
or  phase  angle  has  been  determined  as  before  and  saved. 


the  is  added  to  the  phase  and  the  result  taken  modulo 
N to  determine  the  memory  cell  to  be  accessed. 


As  a simple  example,  let  p (t)  be  uniform  over 


one  sample  period  T . The  transform  of  a pulse  of 

s 


width  Tg  is  the  familiar  Tg  sinc(fTg).  Since  Pd(t) 

0 

is  the  autocorrelation  function  of  p (t),  P.(f)  becomes 

v • ^ 


P^f)  = Tg  sinc2CfTg) 


This  has  its  first  zero  at  f = f » 1/Tg.  Since  the 
power  spectrum  of  the  signal  to  be  resampled  was  shown 
to  be  already  reduced  by  a sine  (fT  ) envelope  due  to 
an  effective  "box  car"  filter,  the  final  spectrum 
before  resampling  has  been  reduced  by  sinc^(fT_)  with- 

5 
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uniform  Jittered  sampling.  More  desirable  filter 
functions  may  be  obtainable  by  careful  design  of  the 
density  function,  but  the  practical  advantage  of  the 
technique  would  soon  disappear. 

In  summary  of  Chapter  4 then,  a technique  for 
generating  harmonically  related  narrowband  processes 
directly  in  the  time  domain  has  been  described. 

The  technique  is  designed  to  be  implemented  in  a digital 
system  where  the  output  must  be  constrained  to  signal 
samples  at  some  system-determined  sample  rate.  The 
spectrum  of  the  resulting  generated  signal  was  described 
and  shown  to  have  potentially  severe  aliasing  problems. 
Increased  memory  for  waveform  storage,  interpolation 
schemes,  and  finally  a statistical  filtering  technique 
based  on  "jittered"  sampling  are  discussed  as  means 
of  reducing  the  aliasing  to  tolerable  levels. 
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CHAPTER  V 

SUMMARY  AND  CONCLUSION 

In  summary,  several  algorithms  for  the  generation 
of  steady-state  signal  components  with  all-digital  systems 
have  been  presented  and  analyzed.  Two  abbreviated  IPFT 
algorithms  for  the  synthesis  of  a broad-band  random  process 
having  a controllable  Power  Spectral  Density  function  were 
analyzed  and  compared  to  a more  conventional  approach  of 
passing  white  noise  through  a Finite  Impulse  Response  (FIR) 
filter.  The  simpler  of  the  two  algorithms  involves  concate- 
nating segments  of  signal  having  an  appropriate  spectral 
composition,  and  results  in  a savings  of  75*  of  the  computa- 
tional effort  required  in  the  FIR  approach.  However,  the 
synthesized  signal  has  a discontinuity  at  each  segment 
boundary  which  may  restrict  its  application.  The  second 
algorithm  employs  summing  an  overlap  of  appropriately 
weight-  -J>  segments  to  eliminate  the  discontinuity,  although 
the  derivatives  at  the  segment  ends  are  still  discontinuous. 
The  overlap  version  results  in  a 50%  savings  in  computation 
load  over  the  FIR  approach. 

For  the  synthesis  of  discrete  components  and  very 
narrow-band  processes,  a pair  of  algorithms  for  generating 
approximation  segments  of  the  desired  signal  by  inverse 
FFT  were  investigated.  Here  the  simplest  version  employing 
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no  overlap  in  output  segments  gives  relatively  poor 
results.  On  the  other  hand,  a corresponding  version 
employing  the  sum  of  appropriately  weighted  and  overlapped 
segments  gives  ve**  good  results  with  a minimal  computa- 
tional effort.  The  two  overlap  algorithms,  one  for 
broadband  components  and  one  for  narrowband  or  discrete 
components,  may  be  combined  to  produce  efficiently 
with  an  FFT  based  system  a signal  having  both  types 
of  components. 

Finally,  an  algorithm  for  producing  harmonically 
rich,  periodic  signals  in  the  time  domain  was  developed 
and  analyzed  and  shown  to  have  potentially  severe 
aliasing  problems.  Various  techniques  to  minimize 
the  undesirable  effects  of  the  aliasing  are  discussed. 

The  result  is  an  efficient,  highly  implementable 
algorithm  for  generating  harmonically  related  signal 
components. 

The  synthesis  of  signals  with  specific  controllable 
¥ 

characteristics  is  a relatively  new  and  still  developing 
part  of  the  signal  processing  field.  As  the  quality 
and  sophistication  of  signal  generation  systems  improves, 
more  and  more  applications  are  found  for  their  outputs. 
Especially  with  all  digital  systems,  the  ability  to 
produce  signals  modeling  problems  with  random  processes, 
and  yet  to  be  able  to  reproduce  the  synthetic  signal 
exactly  or  with  a controlled  perturbation  is  extremely 
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useful  In  evaluating  the  performance  of  a signal 
analysis  system.  The  application  to  simulators, 
particularly  In  training  devices,  is  readily  apparent.  • 
However,  although  many  new  doors  are  opened  by  the 
advent  of  practical,  all-digital  signal  synthesis,  the 

I results  of  this  paper  clearly  indicate  the  need  for 

» 

| careful  analysis  in  applying  the  various  techniques 

| to  any  specific  problem. 

} 

I Several  of  the  results  in  the  paper  have  potential 

i 

! application  in  other  areas  of  digital  signal  processing. 

The  technique  discussed  in  Chapter  3 of  summing  over- 
lapped time  blocks  which  have  been  derived  by  inverse 
transforming  Hanning  weighted  spectra  forms  the  basis  for 
an  extremely  flexible  FFT  filtering  algorithm. 

The  results  found  in  Chapter  4 for  the  spectrum  of 
a resampled  digital  signal  stored  in  a delay  line  or 
memory  is  of  interest,  for  example,  in  a system  where 
a sampled  signal^ passing  through  a shift  register 
delay  line  is  observed  by  a tap,  the  position  of  which 
is  changing  with  time.  Also  in  Chapter  4,  the  concept 
of  using  randomly  disturbed  sampling  times  to  accomplish 
at  least  rudimentary  filtering  as  an  aid  in  minimizing 
aliasing  would  seem  useful,  although  over  a decade 
has  passed  since  Balakrishnan's  paper  and  little  or  no 
application  has  been  made  of  which  the  author  is  aware. 
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An  effort  to  detemlne  what  useful  filter  functions 
night  be  obtainable  under  the  constraints  imposed  on 
its  Fourier  transform  would  probably  be  worthwhile. 
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I APPENDIX  A 

! 


Consider  a truncated  summation  of  the  form 


N-l  „ 

s - l xn 

n*0 


(Al) 


Manipulating, 


N 


XS  * I xn  « XN  + S - 1 
n*l 


(A2) 


S - XS  * 1 - XN 


(A3) 


or 


# 


“ T^X 


(AH) 


Suppose  X is  a complex  exponential  such  that  S is  of  the  form 


S 


N-l 

I 


n*0 


-J^nu 


Applying  AH , with  X « e 


(A5) 
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1 


1 


e 


e 


-j  2*u 


(A6) 


Factoring  an  e“^wU  from 
denominator. 


-jjju 

the  numerator  and  an  e from  the 


(A7) 


(A8) 


Using  the  definition 


8indN(u) 


1 


sin  wu 
sin  Ju 


the  desired  result  is 


(A9) 


N-l  -J  Trtiu 
J e T «Ne  N 
n«0 


sind^Cu) 


(A10) 
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The  Rindjg(u)  function  is  plotted  in  Figure  A1  for  N«9. 
Note  that  for  u equal  to  integral  multiples  of  N,  the 
sind  has  unit  magnitude;  for  all  other  integral  values 
of  u,  the  sind  is  zero.  For  even  N,  the  peaks  alternate 
between  plus  and  minus  one;  for  odd  N all  peaks<  are 
positive  unity.  The  sindN(uf)  is  always  periodic  with 
period  2N;  the  magni-^  ;e  (or  square)  is  periodic  with 
period  N. 


SIND*.  (X) 
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APPENDIX  B 

This  appendix  derives  closed  form  expressions  for 
the  infinite  series 

S,(u,N)  * l sine  (u+rN)  (Bl) 

r*=-- 

and 

• ? 

Sp(u,N)*  l sine  (u+rN)  (B2) 

X»rnm.m 

/• 

that  is,  the  sum  of  a periodic  series  of  sinc(u)  functions 
or  their  squares,  spaced  along  the  u axis  by  a distance  N. 

Consider  first  the  spectrum  G(f)  of  a function  g(t) 
consisting  of  a series  of  N impulses,,  N odd,  with  unit 
spacing  and  centered  at  t»o  as  shown  in  Figure  Bl. 

♦ 

Direct  computation  of  G(f)  gives 

G(f)  « lm  g(t  ) e~J2irft  dt  • (B3) 

rnrnm 

rm  CN-D/2 

G(f ) « / T 6(t-n)  e •J2*rt  dt  (B4) 

-«  n»-(N-l  )/2 

Interchanging  the  order  of  summation  and  integration 


91 


NOLTR  lk-215 


(N-D/2  ■ .... 

G(f)  * I I «(t-n)  e“J2irft  dt 

n*- (N-D/2  — • 


G(f) 


N-l 

TT 


n= 


1) 


e 


-J2»fn 


Writing  out  the  series, 


G(f)  * e d + e ' + ...  + e 


-J2wf^ 


Factoring  out  an  e 


fl(f)  - 1 


N-l  -J2»fn 

" V 


n*o 


Using  A10,  with  u * Nf 


G(f)  - eJ (N-l)  . e-Jwf(N-l)  . N 8lndN(Nf) 


or 


G(f ) « N sind»(Nf ) 
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Now  consider  the  same  function  g(t),  but  derived  as  an 
infinite  Impulse  train 


g*Ct)  ■ l 6(t+n) 

n«-» 


(Bll) 


multiplied  by  a unit  pulse  of  width  N 


p(t) 


i ; 


itiif 


0 ; |t|  » | 


(B12) 


The  spectrum  of  the  impulse  train  is 


a»(f)  * l «(f+r) 


(B13) 


and  the  spectrum  of  the  pulse  is 

P(f)  * N sine (Nf ) (Bl4) 

The  spectrum  of  the  product  is  the  convolution  of 
G»(f)  with  P(f ): 
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OCf)  - /*  G’ (y)P(f-v)dy  (B15) 


0(f)  m r l «(v-lr)  N slnc[(f-v)N]dv  (B16) 

— *»  x^—  ** 


Again  interchanging  the  order  of  summation  and  integration, 

0(f)  « N l J"  sihc[  (f-v,)N]«(v+rOdv  (B17) 


or 


m 

G(f)  « N l sincCfN+rNj 
r*-«» 


(Bl8) 


But  both  BIO  and  Bl8  are  expressions  for  the  spectrum  of 
the  same  function  g(t),  and  thus  must  be  equal.  Therefore, 
letting  Nf  « u, 


l * sinc(u+rN)  ■ S1(u,N)  « sind^Cu) 


(B19) 


To  obtain  the  desired  relation  for  S^(u,N),  consider  the 
auto-conyolution  of  g(t).  By  forming  the  auto-convolution 
of  a series  of  finite  pulses  of  width  6 and  height  1/6, 
and  then  taking  the  limit  as  6 0,  it  can  easily  be  seen  that 

the  auto-convolution  of  g(t)  is  a series  of  2N  - 1 impulses 
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g2(t)  with  a triangular  weighting  as  illustrated  in 
Figure  B2.  The  spectrum  of  g2(t)  is  of  course  Just 
the  square  of  0(f), 


02(f)  - N2sind2(Nf) 


(B20) 


But  g2(t)  can  also  be  derived  by  multiplying  the  impulse 
train  gf(t)  by  a triangular  pulse  p2(t'r) 

N(1  + t/N) ; - N < t <0 

P2(t)  ■ N(1  - t/N);  0.  <t  <N  (B21) 

0 ; jt|  > N 


The  spectrum  of  p2(t)  is 
# 

P2(f)  ■ N2  sinc2(Nf ) (B22) 

since  the  triangular  pulse  is  Just  the  auto  convolution  of 
the  square  pulse  B12,  and  therefore  has  a spectrum  equal  to 
P2(f).  Convolving  P2(f)  with  G»(f), 

02(f)  * /"  P2(f-v)0*(v)dv  (B23) 
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Q2(f ) « r N2  slnc2C(f-v)N]  £ «(v+r)dv 


(B24) 


Interchanging  summation  and  integration* 


02(f)  « N2  £ r sinc2[(f-v)N]«(v+r)dv 


(B25) 


G2(f)  « N2  £ sincZ(Nf+rN) 


Equating  B21  with  B27, 


(B26) 


N2  l 
ra-4 


Binc2(Nf  + rN) 


N2  sind2(Nf ) 


‘ (B27) 


or,  letting  Nf 


£ sinc2(u  + rN)  » S2(u>N)  *sind^(u)  (B28) 
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